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Line Transfer through Clumpy, Large-Scale Outflows: Lya 
Absorption and Halos around Starforming Galaxies 
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ABSTRACT 

We present constrained radiative transfer calculations of Lya photons propagating 
through clumpy, dusty, large scale outflows, and explore whether we can quantitatively 
explain the Lya halos that have been observed around Lyman Break Galaxies. We 
construct phenomenological models of large-scale outflows which consist of cold clumps 
that are in pressure equilibrium with a constant-velocity hot wind. First we consider 
models in which the cold clumps are distributed symmetrically around the galaxy, and 
in which the clumps undergo a continuous acceleration in its 'circumgalactic' medium 
(CGM). We constrain the properties of the cold clumps (radius, velocity, H I column 
density, & number density) by matching the observed Lya absorption strength of the 
CGM in the spectra of background galaxies. We then insert a Lya source in the center 
of this clumpy outflow, which consists of 10^~^ clumps, and compute observable prop- 
erties of the scattered Lya photons. In these models, the scattered radiation forms 
halos that are significantly more concentrated than is observed. In order to simulta- 
neously reproduce the observed Lya absorption line strengths and the Lya halos, we 
require - preferably bipolar - outflows in which the clumps decelerate after their ini- 
tial acceleration. This deceleration is predicted naturally in 'momentum-driven' wind 
models of clumpy outflows. In models that simultaneously fit the absorption and emis- 
sion line data, the predicted linear polarization is ~ 30 — 40% at a surface brightness 
contour of S* = 10~^^ erg s~^ cm~^ arcsec"^. Our work illustrates clearly that Lya 
emission line halos around starforming galaxies provide valuable constraints on the 
cold gas distribution & kinematics in their circumgalactic medium, and that these 
constraints complement those obtained from absorption line studies alone. 

Key words: galaxies: formation - galaxies: absorption lines - galaxies: intergalactic 
medium - ISM: jets & outflows - radiative transfer - scattering 



1 INTRODUCTION 

Deep narrowband imaging has revealed that star forming 
galajcies are surrounded by large (R ~ 100 kpc) Lya halos 
jHavashino et all 12004 ISteidel et all 1201 ll : iMatsuda et"ai] 
|2012| . see also Fvnbo et al. 1999, Ranch et al. 2008). The 
origin of these halos is still disputed, but they likely encode 
valuable information on both the distribution and kinemat- 
ics of cold gas around galaxies. Understanding this so-called 
'circumgalactic' medium (CGM) is vital to our understand- 
ing of g alaxy formation and evolution. 

Zheng et al.l (|2011a ) attribute the presence of extended 
Lya halos to resonant scattering of Lya photons in the 
CGM. In this model, the gravitational potential well of the 
dark matter halo that is hosting the galaxy gives rise to 
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inflows of overdense, ionized, circumgalactic gas. This in- 
falling, overdense gas contains resi dual H I that is opaque to 
Lyg radiation (also see [B arkana fc Loebl 200 3: Santos 2004; 
iDiikstra et al.ll2007l : IHiev et all 120081 : iLaursen et al.ll201l1 ). 
and can scatter a significant fraction of the Lya flux that 
escapes from galaxies into dilTuse halos. 

However, observations of the gas in the circum galactic 
medium of Lyman Break Galaxies (LBGs) indicate that the 
'cold' (i.e. T ~ 10'* — 10^ K) gas is almost exclusively out- 
fiowing: low ionization metal absorption lines are typically 
blue-shifted relative to the galaxies' systematic redshift, and 
the co vering factor of the se blueshifted absorpt ion lines are 
large jSteidel et al.ll2010l '). ISteidel et all l|2010l ') argue that 
these o utfiows -which are n ot present in the the simula- 
tions bv lZheng et al.l (|2011ah - play an important role in the 
Lya transport problem, and scattering through outflows can 
explain (i) in particular the extended red wing of the Lya 
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spectral shape of the Lya hne that is observed in their galax- 
ies, and (ti) the ob s erved Lyg ha l es aro und their galaxies 
l|Steidelet allboill ). ISteidel et all (|2010l ') provide a simple 
model for the Lya scattering process in which a spherically 
symmetric distribution of outflowing clumps surrounds each 
galaxy. The 'covering factor' of clumps on the sky decreases, 
and their outflow velocity increases with distance from the 

galax^^ 

ISteidel et alj (|201Cl ) treat the scattering of Lya photons 
in only an approximate way, and it is not clear whether a 
full radiative transfer calculation would yield similar results. 
Give n the simple geome try of the clumpy outflow model 
that ISteidel et al.l (I2OI0I ) propose, it is straightforward to 
treat the scattering process itself more accurately, and to 
investigate whether their model can indeed quantitatively 
reproduce the observed Lya halos. It is important to test 
the scattering model, and to investigate if other processes 
need to be invoked to explain the presence of extended Lya 
halos around (all) star forming galaxies. These other pro- 
cesses include for ex ample resonant sca ttering in the ion- 
ized inflowing CGM (jZheng et al.ll2011al ). and spatially ex- 
tended Lya emission (a s opposed to sca ttering) from su- 
pernova driven outflows (jMori et al.ll2004l ). or from the cold 
gas streams that provide galaxies with their gas in cosmo- 
logical hydrodynamical simulations (e.g . iFardal et al.ii2001i; 
Furlanetto et al.ll2005l: iDiikstra fc Loeb [20091: iGoerdtet all 



2010l : lFaucher-Giguere et alJl2010l:lRosdahl fc Blaizodl201lD . 
Indeed, recent work has indicated that so-called 'cold-flows' 
may reproduce the observed Lya absorpti on line strengths 
in th e CGM of simulated LBGs quite well l|Fumagalli et al.l 
I2OIII ). which raises the question to whether the Lya halos 
are also related to cold streamy. 

The goal of this paper is simple: to test whether scat- 
tering through a large-scale clumpy, possibly dusty, outflow 
can give rise to spatially extended Lya halos around star 
forming galaxies, and importantly, whether such models are 
consistent with the observed Lya absorption strength of 
the circumgalactic medium as measured in the spectra of 
background galaxies (as in Steidel et al. 2010). In order 
to properly model the scattering process, we modify the 
Mont e-Carlo radiat i ve tra nsfer code for Lya propagation 
from IDiikstra et al.l (|200d). so that it works on arbitrary 
distributions of dusty clumps. 

Having such a code is very useful, as understanding 
Lya transfer through (clumpy, dusty) outflows is relevant 
in a wider range of astrophysical contexts: Firstly, Lya ra- 
diative transfer through outflows affects the imprint that 
reionization leaves on the visibility of the Lya emission line 



^ 'Gravitational heating' of cold flow gas in dark matter halos 
of mass Mjjaio ~ 10^^ Mq, the approximate mass of the dark 
matter halo masses associated with LBGs (e.g. iRudie et al.ll2012l . 
and references therei n), can give rise to Lya luminosities of or- 



der L i^lO"^^ erg s-l JHaiman et al.ll2000l:|piikstra fc Loeljl200g| : 



iGoerdt et aI.ll201Cll : iFaucher-Giguere et aI.l201Cl ). which is an or- 
der of magnitude fainter than the observed luminosities of Lya 
halos around LBGs. This rules out the possibility that gravita- 
tionally heated cold flows contributed signiflcantly to the observed 
luminosity of Lya halos (but see Rosdahl & Blaizot 2011). How- 
ever, it is possible to increase the Lya emissivity of cold stream 
gas if sources embedded in the stream photoionize th e gas in the 
streams (e.g. IDiikstra fc LoeblboogI : iLatif et al.ll20lj) . 



from high redshift galaxies (Santos 2004; Dijkstr a fc Wvithj 
I2OIOI : iDiikstra et al. 201 1; D aval fc Fcr rara 201 1) . Secondly, 
scattering through outflows strongly affects the large scale 
cluster ing of Lya selected galaxies in the post-reio nization 
epoch llZheng et al.ll2011bl : rWvithe fc Diikstrall201ll '). which 
is directly relevant for e.g. the HETDEXJjdark energy exper- 
iment (Hill et al. 20041 . Finally, Lya line transfer through 
clumpy, dusty (not necessarily outflowing) media is of inter- 
est because it can strongly boost the fraction of Lya photons 
that can escape from the interstellar medium of a galaxy, 
possibly even enhancing the equiva lent w i dth of the line as it 
emerges from the galaxy (iNeufelcllll99ll: iHaiman fc Spaang 



I1999I : iHansen fc OH I2OO6I : iFinkelstein et al 1 12009| . but also 
see Scarlata et al. 2009, Kornei et al. 2010). 

This paper is organized as follows: In §[2] we introduce 
the basic quantities that describe a general clump distri- 
bution. In § [3] we describe our model for the cold clumps 
embedded within a hot, large scale outflow, and how we 
constrain the parameters of this model by matching to the 
absorption line data of Steidel et al. (2010). This procedure 
flxes the properties of the scattering medium. We describe 
how we perform Lya radiative transfer calculations through 
this medium in §[4] We present the main results of our calcu- 
lations, and explore how these depend on our assumed initial 
line proflle and dust content of the clumps in § [5] In § 16.11 
we explore more realistic velocity proflles of the cold clumps. 
We discuss model uncertainties and broader implications of 
our work in § [7] before presenting our main conclusions in 
§111 

In this paper we express frequency v in terms of the 
dimensionless variable x = (v — Va)l ti.Va- Here, v^ ~ 2.46 x 
10^^ Hz denotes the frequency corresponding the Lya reso- 
nance (similarly Ac = 1215.67 A corresponds to wavelength 
of this transition), and Av^ = Va\/2kT jm^c?- = lyaVth/c. 
Here, T denotes the temperature of the gas that is scatter- 
ing the Lya radiation. Table [T] gives an overview of sym- 
bols used in this paper. The cosmological parameters used 
throughout our discussi on are (Qm,^A,h) = (0.3, 0.7, 0.7), 
which is consistent with JKomatsu et al.l ([20091). 



2 SOME GENERAL CLUMP STATISTICS 

We denote the number density of clumps at a separation r 
from the source (located at r = 0) by nc{r), and their out- 
flow velocity by Vc{r). We denote the radius of a clump by 
Rc{r). The number density of neutral hydrogen atoms inside 
clumps is denoted by nui{r). The average Hl-column den- 
sity of the cold clumps is denoted by A^hi [r) and is giverjj by 
Nui{r) = 4:nui{r)Rc{r) /3. Finally, the covering factor fc{r) 
is given by fc{r) — nc{r)ac{r), where o-c{r) = TvRc{r)^ de- 
notes the geometric cross-section of a clump of radius Rdj). 
The covering factor fc{r) thus has units of length^ ^, and 
plays a role that is analogous to opacity K(r) in a homoge- 
neous medium. We caution the reader that our deflnition of 



^ http://hetdex.org/ 

^ The average column density is given by Niii(r) = 
^nHr) /cf ^y 27ry7VHi(y), where Nm(y) denotes the total HI 
column density at impact parameter y through the clump, and is 
given by Nuiiv) = 2 yiJ§— j/^njji ■ We can evaluate the integral 
analytically and obtain Niii{r) = 4nHi(f)-Rc(^)/3. 
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Table 1. Summary of symbols used throughout this paper. 



Symbols describing clump properties. 



r 
b 

Tie 
"H 
"HI 
Ra 

rUc 

Nth 






X 

(Ta{x) 

<t>ix) 

dy 

'^in/out 
"in/out 



physical separation of a clump from the source (kpc) 

'impact' parameter (kpe): 

number density of clumps (in kpc~^) 



number density of H-nuclei (i.e. protons, 



') 



number density of neutral hydrogen atoms (in cm~^) 

radius of the clump 

cross-sectional area of the clump (in kpc^): 

(Tc = 7ri?2 

clump mass (in Mq) 

column density of H I through the clump, 

weighted by cross-sectional area (in cm~'^): 

A^Hi = 4iJcnHi/3 

covering factor (in kpc~^): 

/c = nc(Tc 

gas temperature in the cold clumps 

'thermal' velocity of H I in cold clumps (km s~^) 

Vth = ^/2kTc/mp 

Symbols used for Lyct radiative transfer. 

Lya resonance frequency 
Ua = 2.46 X 10^5 Hz 
thermal line broadening (Hz) 

AUa = VaVt\i/c 

dimensionless photon frequency 

X= {v - Va)l l^Va 

Lya absorption cross-section at line center 

o-o = 5.88 X 10-"(Tc/104 K)-i/2 cm^ 

Lycf absorption cross-section at frequency x (cm^) 

daix) = cro4>ix) 

Voigt function at frequency x. 

We adopt the normalization such that 4i{x = 0) = 1 

, and therefore that J <p(x)dx = y/n. 

Voigt parameter 

Ov = 4.7 X 10-''(Tc/10'' K)-i/2 

Unit vector that denotes the propagation direction 

of the photons before/after scattering. 

Unit vector that denotes the electric vector 

of the photons before/after scattering. 

Cosine of the scattering angle 

Scattering phase function 

We adopt the normalization f-^ dfi P(/^) = 47r 



coveri ng factor differs from that adopted by iHansen fc OhI 
l|200Q ). who defined the covering factor to be the mean total 
number of clumps encountered along a random line of sight 
(which we denote by A/'ciump, see below), and which is there- 
fore analogous to optical depth in a homogeneous medium. 
Table [l] summarizes the symbols that we use to describe the 
clump properties. 

Many useful properties of the clumps can expressed in 
terms of this covering factor: 

• The 'mean free path' between clumps Xc{r) — 
l/[nc(r)ae(r)] = 1/Ur). 

• Hence, the mean number of clumps that a random sight- 
line through the distribution of clumps 



where b denotes the 'impact parameter' of the sightline, 
which is the perpendicular distance of the sightline to the 
origin r — 0. We integrate out to radius rmax ~ 250 kpc, 
which corresponds to the radius where the observed absorp- 
tion vanishes (this corresponds to ReS in Steidel et al. 2010). 
• The 'shell covering factor' Fc{r) which denotes the frac- 
tion of the area of a spherical shell at radius r that is embed- 
ded within clumps (denoted by Cf{r) by Martin & Bouche 
2009) is given by 



/■CO 

Fc{r) = / dr' nc{r')acir\r') = 
Jo 

rr+Ra(r) 

■K I dr' nc{r')[R^{r') — {r — r') ] ~ 

Ji — Rc{r) 



(2) 



A/'ciump (6) = 2 



rdr 



&2 



Ur) 



(1) 



g-nc(r)7ri?g(r) = -Ur)R,(r) = /„(0, 

where cjc{r-\r') denotes the area of the clump whose center 
lies at radius r' on the sphere of radius r. The approximation 
on the second line is only true when dRc{r)/dr <^ 1, which 
is generally true in this paper. The quantity fv{r) denotes 
the volume filling factor of clumps at radius r, and is given 
by/„W=nc(r)|^i?c'M. 

If the total mass outflow rate is given by Mout, then 
the total mass fiux through each radial shell is given by 
Mout/'iirr^- The mass densitju in clumps at radius r is then 
Pc{r) — Mout/iTfr^Vcir). If the clumps have a constant mass, 
then the number density of clumps is nc{r) = iVc/47rr^«c(r), 
where A^c is the total rate at which clumps are ejected. For a 
constant Nc, the radial dependence of ac(r) — fc{r)/nc{r) oc 
r^fc{r)vc{r), i.e. Rc{r) oc rfc {r)vc (r). The total number 
of clumps is given by A^ciump = Jq""^'' dr A-Kr'^nc{r). 



3 MODELING CLUMPY OUTFLOWS 
AROUND LEGS 

The theory of large-scale outflows around star forming 
galaxies is extremely complex. Furthermore, the kinemat- 
ics and distribution of cold gas in the outflow is particu- 
larly uncertain. Energy and momentum deposition by ra- 
diation, supernova explosions and cosmic rays create hot 
overpressurized bubbles, which sweep up the surrounding 
IS M into dense, cold shells of ga s (see e.g. the introduction 
of iDaUa Vecchia fc Schavd l2008l : ICeverino fc KlvpinI |2009| . 
and references therein). These 'supershells' accelerate as 
they break out of interstellar medium into the lower density 
CGM, making them subject to hydrodynamical (Rayleigh- 
Taylor) instabilities. The overall acceleration - and hence 
the velocity that the cold gas can reach - depends sensitively 
on when the cold gaseous shells fragment: after fragmenta- 
tion, the hot gas can expand freely through the fragmented 
shell, which reduces the outward pressure on the cold gas 
(e.g. ICooper et al.l l2008l : iFuiita et al.l |2009| . and references 
therein) . 

* The analogy with radiation is illuminating: suppose a central 
source of radiation (instead of mass) has a luminosity L. Then 
the energy flux through radial shell r is s = L/i-Kr^. The energy 
density in the radiation field is u = L/^c-kt'^ (e.g. Rybicki & 
Lightman 1997, also see Fig 1 of Dijkstra & Loeb 2008b). 
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iFuiita et alj (|2009f ) have shown that a spatial resolu- 
tion of ~ 0.1 pc is required to resolve the hydrodynamical 
instabilities, and highlight the physics that likely affects the 
detailed properties of the cold gas in the outflow. For exam- 
ple: magnetic fields may prevent fragmentation of the gas 
shells and thus allow the cold gas to reach larger velocities; 
thermal conduction may further stabilize the cold shells. On 
the other hand, photoionization may heat the cold clumps 
to higher temperatures, which reduces the density contrast 
with the hot wind, which may in turn enhance the fragmen- 
tation of the cold shells. 

Regardless of these complex model details, we expect 
cold fragments of gas to be entrained within a hot wind, 
and that these cold clumps have outflow velocities that are 
less than or equal to the hot win d outflow ve l ocity. The fate 
of these cold clumps is unclear: iKlein et al.l (|l994l ) showed 
that the cold clumps are destroyed on a short time scale as a 
result of hydrodynamical instabilities at the cold cloud-hot 
wind interface. Recent studies have shown that including 
thermal conduction and magnetic flelds in the calculations 
may signific antly enhance the survival probability o f the cold 
clumps (see ICooper et al.l l2008l : iFuiita et al.l |2009| . and ref- 
erences therein). Additionally, new cold clumps might form 
out from thermal instabilities in the hot wind (or hot halo 
gas, e.g. Joung et al. 2012). These two effects (cloud de- 
struction and formation) introduce further uncertainties to 
the spatial distribution of cold clumps. 

The previous illustrates clearly that ab initio modeling 
of cold gas in a large scale galactic outflow is extremely com- 
plex. We therefore adopt a simple pheno menological model 
for th e cold clouds in the outflow as in iMartin fc Bouchg 
l|2009l) . Following Steidel et al. (2010), we assume that these 
clouds are distributed symmetrically around the galaxy. 
This simple model contains parameters, whic h we will con- 
strain by matching the absorption line data of ISteidel et al.l 
(|2010l ). 



individual clump, rric, relates to the total number of clumps 
in our Monte-Carlo simulation as: 



Afc 



dr 



77c SFR 



dr 



(3) 



A'"ciump Jo ^c{r) A'"ciump Jo ^'c(r-) 

where we introduced the cold-gas mass loading factor rjc = 
A/c/SFR. For a given A'ciump we need to know the velocity 
profile Vc{r) to compute the mass in the cold clumps. In 
order to compute the HI gas density inside the clump, we 
need to know press ure in the hot wind. 

Following Mar tin fc Bouchg (12009) we assume a steady- 
state constant velocity hot wind, for which mass conser- 
vation implies Ph{r) = Mh/(47rr^Uh), in which Wh denotes 
the outflow velocity of the hot wind. The number density 
of particles in the hot- wind is nh(r-) « ph{r)/mp, where 
we assumed that the mean molecular weight /.ih ~ 1. Both 
components obey p{r)n '(r) ^constant, where 7 denotes 
the gas' adiabatic index, and their gas pressure is given by 
p{r) = n{r)k-BT{r). Assuming pressure equilibrium between 
the hot and cold gas, we find 

(4) 



Th(r)=Th,o 



rc(r-)=Tc,o 



-27h+27h/7c 



"-c,H 



Tiifi 



Mh 



Tc,o mp 47rr^i^i;h \r.. 



-27h/7c 



Here, 7^ [7h] denotes the adiabatic index of the cold [liotjgas, 
and Tc,o [Tk.o] denotes the temperature of the cold [hot] gas 
at some reference radius r-min, which denotes the 'launch' 
or 'break-out' radius (as in Steidel et al. 2010). We further 
assume that the (constant) velocity of the hot wind is related 
to the temperature of the hot gas at the break out radius a s 
Wh « 940(rh,o/10'' K)'''^ km s~'^ (|Martin fc Bouchell2009l ). 
If we substitute some 'typical' values the we find 



-nl/2 



3.1 The Model 

The total mass outflow rate Mout = r;xSFR, where '77' de- 
notes the 'mass-loadin g' factor (e.g. lOppenheimer fc Davg 
l2006l : lDave et al.ll201lh . We assume SFR=34 Mg yr'^ which 
corresponds to the median UV-based dust-corrected SFR of 
all 92 continuum selected galaxies that were used to create 
the stacked Lya image that revealed the Lya halo surround- 
ing therrlj. We assume the outflow consists of a 'cold' com- 
ponent' which is embedded in a 'hot' component, both of 
which are in pressure equilibrium. The total cold [hot] mass 
outflow rate is denoted by Mc [Mh] . We further assume that 
the clumps all have the same mass (and radius) when they 
are driven out, and that their mass does not change while 
they propagate out. Under this assumption the mass of an 



^ The galaxies that were used to compile the stacked Lya image 
are not the same galaxies for which the CGM was probed with 
background galaxies (see § IT.lj I. However, both sets of galaxies 
were selected in a very similar way and have very similar phys- 
ical properties. For example, the median star formation rate o f 
the 'CGM galaxies' was SFR= 30 Mq yr"! llErb et al.|[2006l) . 
which is practically indistinguishable from the value that we have 
adopted. 



nc,H 



where Th • 



;36 



M, 



\ J- r dl : „ 1 ' ^ / m 1 t-i ' 



r 

"^niin ' 



-27h/7c 



(5) 



Th^o/lO' K, Tc,4 = Tco/lO" K, r^i„,i 



?'min/kpc, and Mh,io = Afh/[10 Mg yr~^]. The number den- 
sity nc,H refers to the total number density of hydrogen 
nuclei in the clump. When the clump self- shields, we ex- 
pect all of the hydrogen to be neutral. Recent hydrody- 
namical simulations of cosmological volumes indicate that 
a decent approximation to the full radiative transfer of ion- 
izing radiation is obtained by assuming that gas self-shields 
when the number density e xceeds some threshold value of 
riciit ^6 X 10~^ cm~^ (e.g. iNagamine et al.ll2010l ). In our 
model, we will assume that the number density of neutral 
hydrogen atoms, nHi,c('") = 7^c,H for nc,H ^ "-crit. When 
?^c,H < n-crit, we assume photoionization equilibrium with 
the UV background, and that the neutral fraction of hy- 
drogen by number is given by xhi = OBn.c,H/rbg. Here, 
QB = 2.6 X 10"^^(rc/10'')"°'' cm^ s"^ denotes the case-B 
recombination coefficient and Fbg = 5 x lO"'^^ s~^ de notes 
the photoionization rate IjFaucher-Giguere et al.ll2008l ). 

To complete the description of our outfiow model, we 
need to assume Vc{r). As discussed previously, this velocity 
profile is not well known, and depends sensitively on when 
the cold gas shells fragment, and on whether cold clouds 
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form from the hot wind as a result of thermal instabilities. 
We follow the empirical approach of Steidel et al. (2010), 
who assumed that the acceleration of the cold clumps scales 
as ac{r) — Ar~" , which results in a velocity profile of the 
form 



«cM=(^)"^r--.-)° 



(6) 



for a > 1 (|Steidel et alj |2010|), where A is a constant 
that sets the velocity at r — > cxd, v^d- That is, Hoo = 

y'^Ar^^^ /{a — 1). We do not consider models with q ^ 1 
in this paper (see belowjj. 

Formally, our model has thus ten parameters which in- 
clude: rmin, a, Voa, Th.o, Tc,o, 7h, 7c, ??c, ??h and iVciump- For 
most of these parameters we have decent constra ints, and 
they are not free. For example, ISteidel et al.l (|201Cl ) inferred 
from their observations that 1.15 < a < 1.95, i)cx> = 800 km 
s~^ and adopted rmin = 1 kpc, w hich is close to theoretica l 
estimates of the blow-out radius (jMartin fc Bouchg 120091 ). 
The observationally inferred hot wind temp eratures lie i n 
the range Th = 10*^ - lO'^ K for dwarf galaxies (fMartin"l999"), 
but could be larger by a factor of ~ 10 in starburst galaxies 
fe.g. lStrickland fc Heckmanl2009h . Cold gas at temperatures 
Tc > 10* K would efficiently cool down to Tc ~ 10* K, below 
which gas cooling becomes less efficient. Further cooling is 
possible because of metals, but it is unclear to which tem- 
peratures the gas can cool in the cold clumps, and to what 
extent the clumps are heated as a result of their interaction 
with the hot wind (and/or as a resul t of photoheating) . We 
will co nsider values of log Tc = 2 — 4. lOppenheimer fc Davg 
(|2006) could reproduce the observed amount of C IV absorp- 
tion in quasar spectra at z=2-5 with large scale (momentum- 
driven) outflows arising from star forming galaxies, for a 
total mass loading factor ry = r/c -I- ?7h = aoja. Here, g de- 
notes the velocity dispersion of the galaxy, and uo = 150 km 
s~^ (|Qppenheimer fc Davgl2008l ). The observed dispersion 
of the Ha line in the galaxies that were used to construct the 
stacked Lya image, is cthcv ~ 100 — 150 km s~^ (see Fig 4 of 
Steidel et al. 2010). Under the reasonable assumption that 
ana provides a good measure of a, we require that the total 
mass loading factor is close to unity. However, we caution 
that direct observational constraints on r\ a re uncertain by 
a factor of at least a few (JGenel et al.ll2012l . and references 
therein). The adiabatic index of the cold gas is 1 ^ 7 ^ 5/3, 
where 7 = 1 [7 = 5/3] corresponds to isothermal [adiabatic] 
gas. 

In this paper, we only consider isothermal cold clouds, 
i.e. 7c = 1, because we found that for models with 7c > 1.2 
generally the cold clumps are compressed too much, which 
gives rise to EW-6 curves that decline too steeply with h. 
We also only consider models with either 7h = 1 and 7h = 
1.2. As we will show, in models with 7h = 1.2 the cold 
clumps expand more rapidly, which causes nc,H < ^crit, and 



the HI column density of the cold clumps declines too fast 
for significant scattering of Lya photons. Finally, we only 
consider A'ciump = 10^ and A'ciump = 10^. This choice for 
A'ciump translates to clump masses in the range ttIc = 10* — 
10^ M0 (see Table [2|. We stress that our main results are 
insensitive to the adopted value for A'ciump- 

3.2 Constraining the Free Parameters of the 
Model with Absorption Line Data 

As was discussed above, each outflow model contains ten 
parameters, seven of which are allowed to vary within 
a reasonably well known range. Each model is there- 
fore parameterized by the 7-D parameter vector P = 
('■min,a,Uoo,T'h,o,rc,o,?7c,??h). For each P we obtain a dis- 
trib ution of cold c l umps that contain neutral atomic hydro- 
gen. ISteidel et al.l (|20ld ) have measured the average Lya 
absorption line strength at impact parameter h from galax- 
ies by analyzing the spectra of background galaxies. We can 
use these observ ations to cons t rain t he components of P. 

Specifically, ISteidel et al.l l|2010f ) have measured mean 
absorption equivalent width (rest frame) in the Lya line at 4 
impact parameters. These measurements are shown as filleA 
blue circles in Figure [T] The equivalent width (at impact 
parameter b) is defined as 



EW(6) = A, 



Afa 



dx(l-exp[-T(x,b)]). (7) 



The integral is over the dimensionless frequency x (see Ta- 
ble [l|, and exp[— r(x, b)] denotes the fraction of the flux at 
frequency x that is transmitted to the observer. 

Evaluating this transmission is more complicated for 
a clumpy medium than for a homogeneous medium. For 
example, in the hypothetical case of A/'ciump(fe) =0.1 we 
expect only 10% of all sightlines with impact parame- 
ter 6 to pass through a clump, and the outflow is trans- 
parent for the remaining 90% of the sightlines. The ab- 
sorption line strength is then EW(6) =EWciump('")/10, 
where EWciump('') denotes the equivalent width as a re- 
sult of absorption by a single clump. In the more gen- 
eral case, the transmission exp[— T(j:,fe)] is a product of 
the transmission by individual segments along the sight- 
line. That is exp[— t(2;,6)] = J^^J^ exp[— r(3;, 6, Si)]. Here, 
exp[— r(a;, b, Si)] denotes the fraction of the flux that is trans- 
mitted by clumps in line segment 'i', that lies at line-of- 
sight coordinate si, which has length Asi, and which lies a 
distance ri = ^/b^ + s^ from the galaxy. This transmission 
exp[—T{x,b,Si)] = Pciump,iexp[— rciump(a;')] + 1 — Pdump.i- 
Here, Pciump denotes the probability that a clump lies on 
line segment 'i', and exp[— rciump(a;')l denotes the total frac- 
tion of the flux that is transmitted by the clump. Note that 
because of the outflow velocity of the clump, x' is related to 
X via a Doppler boost (see below). Theprobability that line 
segment 'i' contains a clump is given bjOpciumpa = fc{r)Asi, 



^ In the simulations of iDalla Vecchia fc Schayg 1 12008^ . the wind 
velocity increases with r simply because the gas at a given radius 
has a velocity that is close to the minimum velocity it must have 
had, to reach that radius in the finite time since the launch of 
the wind. Having Vc{r) increase with r does therefore not solely 
represent models in which the clumps accelerate with radius. 



^ Formally, we do not allow the clumps to overlap in our 
Lya Monte-Carlo calculations, which modifies the probability 
Pciump, i = fc{r)Asi X p(no clump at As ^ 2_Rc) ~ /c(r)Asi X 
[1 — 4/c(r)i?c('')]. We have explicitly checked that this difference 
makes no difference in practice because of the low volume filling 
factor for the cold clumps. 
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and we obtain 



JVa 



10 



exp[— r(x, b)] = 



Yl (^As.Mn) exphiVc,Hi(^Oo-«(^')] + 1 - As,/c(rO) . (8) 
1=1 

Here, line segment 'i' covers the range [smin + (i — 

l)ASi,Smin + iASi] , where Smin = — V^max ~ b'^ Billd ASi = 

2|smin|/Afs- Furthermore, aa{x') denotes the Lya absorp- 
tion cross-section evaluated at frequency x' in the frame of 
the clump, i.e. x' — x — ""^''' bb • er. Here, Gs [erl denotes a 
unit vector along the line of sight toward the galaxy [along 
the radial vector]. We stress that the background galaxies 
do not provide sight-Zines through the CGM of the fore- 
ground galaxy. Because of their finite physical sizes, they 
instead probe 'sight-cylinders' of some radius rcyi through 
the CGM. However, it is easy to show that Eq[H]also applies 
to cylinders of radius rcyi provided that rcyi <C 6, which is 
the case for t he observations a t b >30 kpc since rcyi is ~ a 
few kpc (e.g. iLaw et al.ll2012l ). Formally, this formalism is 
not correct at 6 ~ 0. However, we have explicitly verified 
that a more detailed calculatior|f| reduces the predicted EW 
by only ;^20% at 6 = for rcyi ^ 5 kpc. 

To find models that fit the data of Steidel et al. (2010) 
we use a Markov-Chain Monte-Carlo simulation to probe 
the parameter space spanned by P. Our exploration of the 
parameter space is rather simple: our goal is to find models 
that provide a good fit to the data, whether these models 
are physically plausible in the context of our model, and 
to explore whether these same models can give rise to the 
observed Lya halos. We will not present probability distri- 
bution functions for the elements of P, and will not explore 
the correlations that exist between them. Given the simpli- 
fied nature of the model, this would distract from the main 
purpose of our analysis. 

From the Markov-chainqj we select three models, which 
we denote with model I-model III. We summarize the 
parameters of these models in Table [21 Figure [T] compares 
the observed EW as a function of 6 to the predictions by 
the models. The black solid line, red dashed line, and blue 
dotted line represent model I, model II, and model III 
respectively. All models clearly provides a good fit to the 
data. 

Figure [2] shows some properties of the cold clumps for 
all three models, where we use the same line style and color 



* In the case of a cylinder, we replace Pclump.i with the frac- 
tion of the area of the background galaxy, A^g, that is cov- 
ered by clumps in cylinder segment 'i', which is given by 

^ "' /J""*"' dx X n^{u)Ac{u), where u^ = x^ -\- s?. 
^ We generate 10 chains that contain 2000 steps and simply se- 
lect the best-fit model from all ten chains. For each step we 
compute the likelihood C{F) = exp(— x^/2)Pprior(P)i where 
X^ = Ei=l(EWmod,i - EWobs,i)^/o-|w,i- We assume flat priors 
on Tc, rjij, r}c, but restrict ourselves to the range 2 < logTc,o < 4, 
0.1 < T^h < 10, 0.1 < rjc < 10 (i.e. outside this range, the 
prior probability is set to zero). We assume Gaussian priors for 
rmin [(r,fT) = (1.0,1.0)], a [{a,a) = (1.5,0.4)], v^ [{v,a) = 
(800.0,100.0)], and logTh,o [(T, tr) = (7.0,0. 5)], but restrict our - 
selves to < r^in < 3.0 kpc, 0.5 < a < 2.0 llSteidel et al.ll2010t) . 
500 <Vaa < 1000 km s"!, and 5.0 < logTh,o < 8.0. 






m 



0.1 



0.01 







r 


• Steidel-hlO 

model I (text) 

model 11 

model 111 


- 




1 1 1 


" 



10 100 

Impact parameter b (kpc) 



Figure 1. This figure shows the mean absorption line strength - 
quantified by the restframe equivalent width- in the Lya line as 
a function of impact parameter b. Blue filled circles indicate the 
observations by Steidel et al. (2010). The three lines indicate EW 
as a function of b for the three models (model I-III) for the large 
scale outflow (see text for details on the model) . These models are 
used as input to our Lya radiative transfer calculations. 



representation as in Figure [T] We discuss the clump proper- 
ties in more detail below: 

• Model I: The clump radii increase from Re < 0.01 kpc 
to Re ~ 0.1 kpc at r=100 kpc. The HI column density of 
the clumps falls from A^cHI ~ 2 x 10^'' cm"^ at r=10 kpc to 
A^Hi ~ 10^^ cm~^ at r=100 kpc. The lower left panel shows 
that the HI number density stays above ricrit all the way 
out to r ~ 200 kpc, and that then the HI number density 
drops off fast. The central lower panel shows that the outfiow 
velocity of the clumps increases rapidly at small radii, and 
that it barely increases further at r J>10 kpc (this plot is 
also shown in Fig 23 of Steidel et al. 2010). Finally, the 
lower right panel shows that a random sightline at impact 
parameter b intersects ~ 1 clump out to 6 = 100 kpc, after 
which it decreases rapidly. 

• Model II: Most differences between this model and 
model I are easily understood: the number density of 
clumps is lower by a factor of 10 as a result of the lower to- 
tal number of clumps. In order to yield the same absorption 
line strength, the decrease in nc{r) is compensated for by an 
increase of their radii {upper central panel). The HI number 
density in the cold clumps is set entirely by the properties in 
the hot gas, and hence remains identical {lower left panel). 
As a result of the unchanged HI number density and the 
enhanced clump radii, the total HI column density is corre- 
spondingly larger {upper right panel). Finally, because of the 
enhanced HI column density of individual clumps, A/'ciump(6) 
must be smaller for model II in order to reproduce the ob- 
served EW(6) {lower right panel). 

• Model III: The number density of clumps is the al- 
most identical to that of model II, which is because this 
number density is determined mostly by the total mass out- 
flow rate of cold gas (comparable for both models) as well 
as the velocity profile (also comparable for both models). 
The temperature of the hot gas decreases with radius from 
Th ~ 10^ K to Th ~ 10^ K in this model {upper right panel), 
pressure equilibrium also implies that the number density 
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Table 2. Outflow model parameters that we adopt in our models. 
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Figure 2. The plots give a more detailed look at models I-III of the cold clumps in the large scale outflow. We show the number 
density of clumps nc(r) in proper kpc"'' in the upper left panel. The upper central panel shows the clump radius in proper kpc as a 
function of r. The upper right panel shows the HI column density (in cm~^) through an individual clump located at radius r. The lower 
right panel shows the number density of HI atoms in the clump (in cm~^), while the lower central panel shows the velocity profile (in 
km s~^). Finally, the lower right panel shows the number of clumps a random sightline with impact parameter b intersects. 



of HI atoms decreases faster in the model (lower left panel). 
Because we fixed the clump mass, the clump radius is there- 
fore increasing more rapidly with radius in model III than 
in model II (upper central panel). 

3.3 Generating a Random Realization of Clumps 

We generate random realizations of model I-model III 

Eis follows. For each of the A^ciump clumps in our model, we 
first generate a random unit vector ei. Then we generate a 
random radial coordinate ri of the clump from 



7^: 



/' 



dr 4Tvr nc(r), 



(9) 



where 7?. is a random number between and 1. The center 
of the clump 'i' is then given by Xi = nei. We then check 
whether clump 'i' overlaps with any of the previously gener- 
ated i-1 clumps. In case it does, we generate a new random 



unit vector ei, until clump 'i' does not overlap with any of 
the other clumps. We then proceed to generating the posi- 
tion of clump 'i-|-l'. The velocity, temperature, radius, HI 
number density, and dust content are specified fully for a 
given r-i and a given model (the dust content will be dis- 
cussed in 5 15.21) . 



Once we have random realizations for each model, we 
shoot random sightlines at a range of impact parameters 
(100 sightlines at each impact parameter) and compute the 
mean Lya absorption strength as a function of impact pa- 
rameter. We show the results of this exercise in Figure [S] 
where we use the same line colors and style to represent 
the different models as in the previous plots. Our discrete 
realizations of clumps also give rise to Lya absorption at 
levels that are in excellent agreement with the data. This 
provides an independent check of the accuracy of Eq[S]and 
that we have generated representative random realizations 
of our best-fit models. 
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Figure 3. Same as Fig[T] but for random realizations of model 
I— model III. This Figure illustrates that we have generated rep- 
resentative random realizations of our best-fit models. It also pro- 
vides an independent check of our predicted EW as a function of 
b in Eq|8] 



4 THE LYq transfer MODEL 

4.1 Monte-Carlo Calculations 

Our Monte-Carlo Lya radiative transf er code was developed , 
tested and described in more detail in lDiikstra et al.l l|200(i ). 
For our current work we modify the code in several ways: 

• We can follow the transfer through an arbitrary dis- 
tribution of spherical clumps. Each clump is assigned a lo- 
cation, a radius, a hydrogen number density, a dust num- 
ber density, a temperature, and a (ra dial) velocity. This di f- 
fers from the applications presented in lDiikstra et al.l l|2006l ). 
w here we focused on sing le spherically symmetric gas clouds. 
In iDiikstra et al.l l|2006l ) the single clump was divided into 
a large number of spherically concentric shells, to each of 
which we assigned a velocity, HI density, and temperature. 
In this work, we do not allow for gradients of temperature 
etc. across the clump (although our code can be modified 
for this purpose). 

• We include the impact of dust on the radiative trans- 
fer process following Laursen et al. (2009). We assume 
an 'SMC type frequency dependence of the dust absorp- 
tion cross-section. However, this frequency dependence is 
practically irrelevant over the narrow range of frequencies 
which are covered by Lya photons (see Laursen et al. 2009 
for details). When Lya photons scatter off a dust grain, 
we assume that the scattering is described by a Henyey- 
Gree nstein phase function with asymmetry parameter g = 
0.73 iLaursen et al.ll2009l . and references therein). Scatter- 
ing of UV radiation by SMC t ype dust give s rise to little 
linear polarization (see Fig 5 of iDraind 120031 ) . and for sim- 
plicity we shall assume that scattering by dust grains does 
not polarize Lya radiation. For the work presented in this 
paper we will assume that the albedo, which denotes the 
ratio of the scattering to the total cross-section, is A = 0.0 
(see § 15. 2p . We assume that the relation between dust ab- 
sorption opacity ro.a and th e color excess Eb-v i s given by 
T"D,a = lO.O-Bs-v (also see iVerhamme et al.ll2006l ). 

• Because we want the option to study clumpy out- 



flows that are not spherically symmetric, we generate 
surface brightness proflles with the so-called 'peeling 
algorithm', or the 'next event estimator'. This technique 
has been employed in many previous stu d ies of Lya trans- 
fer (e.g^ Zheng fc Miralda-Escudel 20021: Cantalupo et al.l 



"Ml 



20051: iTasitsiomil |2006| : 



Faucher-Ciguere et al 



Laursen fc Sommer-LarsenI 
2010l : 



Kollmeier et al 



2007 : 



20ia : 



Barnes et al.l [201l|). A more detailed description of 



how we ge nerate images can be found in Appendix lAl. 21 

• As in lDiikstra fc Loebl (|2008) we compute polarization 
of scattered radiation. Polarization calculations have thus 
far focused solely on spherically symmetric gas distributions. 

We point out that our code allows us to perfectly 're- 
solve' the Lya transfer inside of clumps with arbitrary small 
radii (in this case Re <\Q pc) (see § 13. 2p within our large 
(diameter ~ 500 kpc) outflow. 



4.2 Analytic Calculations 

Under the assumption that the Lya photons scatter only 
once-which is reasonable as we will argue below- we can 
compute the surface brightness profile [5(6)] as well as the 
polarization profile ['P(6)] analytically as 



V{h) 



S{h) = Si{h) 
Si{b)- 



5.(6)1 



Sl{h)^Sr{h) 



, where Si(h) and 
(JRvbicki fc Loebl 1 19991 : 
are given by 



Sr (b) denote p o larize d 
iDiikstra fc Loebl l2008l ). 



Siib) 

Srib) 



/ kpc 
Vasec 



dx 



as - X - X 
, b 4 



(10) 

(11) 

fluxes 
which 



(12) 



Sin(a;,fe, s) X fc{s,b) x (l - exp[-rciump(a::', s,6)]) x 



xfascix",b,s) X 



where the term (kpc/asec) 
into units erg s~^ —-2 



converts the surface brightness 
cm " arcsec"^. Furthermore, Sin{x,b,s) 
denotes the total incoming flux at location (fe, s) and fre- 
quency X, where s denotes the line-of-sight coordinate (see 
Fig lBip . We can write Sin{x,b,s) = soj^exp{—r[s,b,x]), 
in which so denotes the total observed Lya flux of the 
source if no scattering occurred at alo, r = \/s^ + b^, and 
exp(— r[s,6, x]) denotes the total fraction of the flux at fre- 
quency X that has not been scattered out of the line of sight 
yet. The term exp(— r[s,6, x]) is computed as in Eq[8l The 
optical depth Tciump(a;', s, b) = Niii,c{r)o-a{x'), denotes the 
optical depth through the clump at radius r at frequency x. 
In the frame of the clump, photons of frequency x appear at 
x' = x—Vc{r)/vth- The last line contains the scattering phase 
function, in which n — —s/r, and accounts for the fact that 
photons are not scattered in all directions with equal prob- 
ability. Finally, f csc{x" ,b,s) is the probability that photons 
that are scattered towards the observer, are detected. This 
probability can again be computed as in Eq|51 but note that 

^^ The flux so relates to the intrinsic luminosity. La, of the source 
simple as sq = La47rd£(2), where di^{z) denotes the luminosity 
distance to redshift z. 
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after scattering the photon's frequency has been changed to 
x" = X + {jj. — l)vc{r)/vth- We present a complete derivation 
of this equation in Appendix [Bl 



5 RESULTS 

In our Monte-Carlo calculations we emit A^^ = 10"" Lya pho- 
tons in the center of the clumpy outflow for each model, and 
randomly draw the initial frequency of each emitted photon 
from a Gaussian with a standard deviation of a = 150 km 
s~^. This is close to the typical dispersion of the observed 
Ha lines in the sample, which is about aua ~ 100 km s~^ 
(see Fig 4 of Steidel et al. 2010). We assume that the lumi- 
nosity of the central source is La = 3.4 x 10*'^ erg s~^. This 
luminosity corresponds to the intrinsic Lya luminosity of a 
gala:xy that is forming stars at a rate SFR=34 Mq yr~^, 
which corresponds to the median UV-based dust-corrected 
SFR of all 92 continuum selected galaxies that were used 
to create the stacked Lya image. We thus implicitly assume 
that the escape fraction of Lya photons from the dusty inter- 
stellar medium into the large scale outflow i^^ fcac = 100%. 
The predicted surface brightness scales linearly with /esc. 
The luminosity La = 3.4 x 10''^ erg s~-^ which at z = 2.65 



translates to so = 5.9 x 10~ erg s~ cm~ , which is relevant 
for our analytic calculations. 

Throughout, we represent the observed Lya surface 
brightness profile by the function S{b) = So exp[— 6/bc], 
where b denotes the impact parameter in kpc. Here, So = 
2.4 X 10"^* erg s"^ cm"^ arcsec"^, and 6c = 25.2 kpc. This 
function provides an accurate fit to the average Lya ob- 
served in the full sample of 92 continuum-selected galaxies 
with Lya imaging (Steidel et al. 2011 ), and is shown as red 
solid lines in the Figures below. 

We represent results from our Monte-Carlo calculations 
with black filled solid circles, which contain errorbars. We 
obtain these points by averaging the six surface brightness 
(and polarization) profiles that we obtain by viewing the 
clump distribution from six different directions (see §|47l|. 
Uncertainties on the these points indicate the standard de- 
viation from this average. 



5.1 Dust Free Clumpy Outflows 

The top panels of Figure|3]compare the observe d Lya surface 
bright ness profile (in erg s~^ cm~'^ arcsec"^) of lSteidel et al.l 
l|2011l ) with the 'predicted' surface brightness profiles ob- 
tained from the Monte-Carlo {filled black circles) and ana- 
lytic {green solid Imes) for model I-model III. These fig- 
ures show clearly that none of the models can reproduce the 
observations: all three models have too much flux coming 
from d < 2 arcsec, although the actual observations also 



" The conversion La = W^^ x [SFR/(iWQ/yr)] applies for a 
Salpeter stellar initial mass function (IMF) and solar motallicities. 
We expect a larger Lvq lum inosity at fixed SFR at lower gas 
metallicities l lSchaererl 120031) . For a Chabrier IMF - and again 
solar metallic ity - we expect ~ 1 .8 times more Lyo luminosity at 
a given SFR llSteidel et al.ll201Cr) . It may therefore be possible to 
get the same Lyo luminosity from the central source if the escape 
fraction is /esc ~ 50%, which is still large. 



show a significant excess over the exponential fitting func- 
tion at these inrpact parameters. This difference can be re- 
duced by including dust (see § 15. 2p . The most important 
difference however, is at large impact parameters {6 ^5 arc- 
sec, or b ^40 kpc), where our predicted surface brightness 
profiles are low by an order of magnitude. 

The fact that our model surface brightness profiles are 
so much fainter is easy to understand: the central lower panel 
of Figure [5] shows that at r > 10 kpc, the cold clouds are 
moving away from the central Lya source at Vc <;600 km s~^ . 
The majority (95%) of Lya photons will therefore enter the 
clumps with an apparent redshift of Av = 600 ±2a — 300 — 
900 km s~^. In order for the clouds to be opaque (r > 1) to 
Lya photons, we must have A''hi ~ 3 x lO'^® — 3 x lO^'^ cm~^. 
For most of our models, this requirement is met. However, 
the observed Lya absorption line strength at b = 30 kpc 
(see e.g. Fig[T} is EWobs(6 = 30 kpc) ~ 2 A, which is smaller 
than than the absorption equivalent that would be produced 
by a single clump with this HI column density. The number 
of such clumps along the line of sight is therefore small (order 
unity, see the lower right panel of Fig [5} , and because the 
radii of the clumps Re <C r, the majority of photons escape 
from the circum galactic medium without encountering the 
outflowing, cold clumps. 

This last point is also illustrated nicely by the fraction 
of photons that never scatter in the outfiow, denoted by /ns: 
in model I we find that /ns ~ 0.75. That is, the majority of 
photons do not scatter off the cold clumps that give rise to 
the observed absorption. For model II, we have /ns ~ 0.80, 
while for model III we have /ns = 0.43. An immediate 
implication of this finding is that the photons that do not 
scatter should be observed as a Lya point source of equal or 
larger luminosity than that of the Lya halo. In the obser- 
vations, the luminosity of the Lya halo exceeds that of the 
central source by about a factor of ~ 5, also in disagreement 
with our model. 

Those photons that do scatter in the outfiow, get 
Doppler boosted to lower frequencies after they escape from 
the clump. The probability that the photon scatters in a sec- 
ond clump is then reduced further, because they appear even 
further from line resonance in the frame of the other clumps. 
Indeed, we find in the Monte-Carlo simulations that photons 
generally scatter only in one clump, and this is the reason 
why our analytic solutions for the surface brightness profiles 
closely matches the ones we obtained with the Monte-Carlo 
methoco. The lower panels shows that the scattered Lya 
radiation is highly polarized, with the linear polarization 
V ^40% at 5" ^10"^** erg s"^ cm"^ arcsec"^ This high 
level of polarization is another consequence of the fact that 
the photons typically scatter only oncqij. 



^■^ We stress that we do not expect perfect agreement, mostly 
because: («) in our dust-free Monte-Carlo calculations the total 
flux of Lya photons through each radial shell is conserved (but 
redistributed along the frequency axis), while this is not the case 
for the analytic calculations; {ii) the analytic formulation does 
not properly account for radiative transfer effects inside the clump 
where the photon scatters; {Hi) we construct images that contain 
pixels that are 5x5 kpc^ on a side, and binning is involved when 
creating azimuthally averaged surface brightness profiles. 
^'^ The HI column density of clumps A^hi <;10^^ cm~^ at r <;10 
kpc. Photons typically scatter several times on the surface of these 
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Figure 4. This Figure shows a comparison between the 'predicted' and observed Lya line profiles for our three models. The red solid 
lines show the surface brightness profile that was observed by Steidel et al. (2011). The solid green lines show the surface brightness 
profiles as given by Eq 12 which assumes that photons scatter only once. The filled black circles show the results from our Monte-Carlo 
calculations. This Figure shows that the scattered Lya radiation in our models gives rise to Lyo halos that are more compact (i.e. 
centrally concentrated) and fainter by ~ an order of magnitude at 6 >5 arcsec. The lower panel shows the linear polarization as a 
function of impact parameter. We find that scattering through clumpy outfiows can give rise to high levels of polarization. This figure 
also shows that there is very good agreement between our analytic and Monte-Carlo calculations, which is a consequence of the fact that 
photons typically scatter in only one clump (or none at all). 



Figure [S] shows how our predicted surface brightness 
profiles change when we increase the width of the Lya spec- 
tral line of the central source by a factor of three to cr = 450 
km s~^. Increeising the line width enhances the fraction of 
photons that enter the outflow with a significant blueshift. 
These photons appear closer to the line resonance in the 
frame of the outflowing gas, and are therefore more likely 
scattered. These plots show that broadening the initial line 
proflle increases the surface brightness level of the scattered 
Lya halos, but not nearly enough to account for the ob- 
served Lya halo proflles. Choosing even broader lines would 
clearly increase the surface brightness of the Lya halos more, 
but this choice would probably be unphysical. Scattering of 
Lya photons through large columns of interstellar HI gas 
(A^Hi ^lO^'^ cm~'^) could in theory easily broaden the line 
even more than this, but the line broadening as a result 
of scattering is limited by dust in the ISM (and it's dis- 
tribution, see Figure 8 of Laursen et al. 2009). Moreover, 



clumps before escaping from them. This suppressed our the po- 
larization signal that we predict with the Monte-Carlo calcula- 
tions at small impact parameters. We do not completely 'wash 
out' the polarization signature because polarization measures the 
anisotropy in the local Lya radiation field weighted by the pho- 
tons' escape probabilities. This can lead to significant levels of 
polarization d espite the fact that ph otons can scatter multiple 
times (also see lDiikstra fc Loebll2008l V 



H I column densities of this magnitude would efficiently 
trap Lya radiation. The radiation pressure exerted by this 
trapped radiation can c ause the HI gas to expand outwards 
(JDiikstra fc Loebll2008l ). Scattering of Lya photons by this 
outfiowing gas woul d broaden, but also r edshift the line (e.g. 
lAhn fc Lee 2002; V erhamme et al.ll2006l) . Especially for H I 
column density of A^hi ~ lO^'^ cm~'^, this redshift can be 
large even for outflow velocities of only a few tens of km 
s~^. These redshifted Lya photons would appear further 
from resonance in the frame of the cold clumps, and they 
would less likely be scattered. We therefore think that our 
choice a — 450 km s~^ corresponds to a reasonable upper 
limit on the amount of line broadening that occurs in the in- 
terstellar medium. Finally, the filled black squares show that 
the predicted surface profiles are affected most strongly by 
dust in the central regions (see § I5.2[l . 



The fundamental reason why we fail to reproduce the 
observed Lya halos is that the clumps at large radii are 
moving away from the Lya source too fast, thus requiring 
prohibitively large HI column densities in order to remain 
opaque to the Lya photons. In § 16.11 we explore a class of 
models for which the cold clumps decelerate after some ra- 
dius, as expected naturally in models of momentum driven 
winds. 
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Figure 5. This Figure shows the same as Figure [4] Here, we include models in which we assumed the initial Lya line profile to be 
broader by a factor of 3, i.e. a = 450 km s~^ . The resulting surface brightness profiles and polarization profiles are shown as the black 
dashed lines. We have shown results from our analytic calculations only. These plots show that significant (and probably unrealistic, see 
text) broadening the initial line profile increases the surface brightness level of the scattered Lya halos, but not nearly enough to account 
for the observed Lya halo profiles. The filled black squares show the result of models in which we include (a generous amount of) dust. 
Dust affects the inner regions of the surface brightness profiles most strongly. 



5.2 Dusty Clumpy Outflows 

In all models, the central clumps have the largest column 
densities. If the dust opacity of the clumps scales with 
their HI column density, then we can suppress the ob- 
served flux from the central regions. We can therefore 'flat- 
ten' the predicted surface brightness profile by adding dust 
to the clumps. We investigate this in more detail here. 
We rerun model I-III but add dust to clumps. The to- 
tal amount of dust in a clump is normalized such that 
average dust absorption optical depth through a clump is 
Tdust.aba = fcdust X (A^Hi/10^° cm~^). That is, for a column 
density of A^hi = 10^" cm~^ , the total dust absorption opac- 
ity is Td = fcdust- We will explore the impact of dust for 
fedust = 1. For comparis on, the diffuse HI pha se of the Milky- 
Way has fcdust ~ 0.1 IjHansen fc OhI |2006| . and references 
therein). We choose the larger value fcdust = 1 to clearly 
illustrate the potential impact of dust. To maximize the im- 
pact of dust, we also assumed that the dust grains do not 
scatter Lya photons (i.e. A = Q). 

The black filled squares in Figure [5] show the predicted 
surface brightness and polarization profiles in the presence 
of dust. We indeed find that the surface brightness profiles 
are suppressed most at small impact parameters. Also, the 
predicted polarization is affected very little. Note that for 
non-zero dust scattering albedo, the radiation that was scat- 
tered by dust would be unpolarized, which would lower the 
overall polarization. However, if the Lya halos were a re- 
sult of scattering by dust grains, then we would expect the 



UV-continuum to follo w the same surface brightness profile, 
which is not observed (jSteidel et al.ll201ll ). 



6 LYa HALOS FOR MORE 'REALISTIC 
MODELS 

6.1 Model IV: Adopting a More Realistic 
Velocity Profile 

In our previous models, the cold clouds accelerate as they 
break out of the interstellar medium. The acceleration de- 
creased with radius as adr) oc r~" , where a = 1.4 — 1.5 
in model I-III. This continuous acceleration represents a 
'momentum-driven' wind scenario in which the cloud's ac- 
celeration is driven by for example ram pressure of the hot 
wind or radiation pressure. However, galaxies populate the 
centers of gravitational wells, and the cold clumps are sub- 
ject to a gravitational force which decreases as ex r~^ in 
the potential of an isothermal sphere. The deceleration of 
clumps as a result of gravity therefore decreases slower with 
radius than ac{r), and gravity dominates beyond some some 
'transition' radius rtrans. This deceleration following accel- 
eration occurs for any model in which a > 1. This decelera- 
tion potentially has important implications for the predicted 
surface brightness profiles of the Lya halos: one of the main 
reasons model I-III significantly underpredict the Lya sur- 
face brightness is that the clumps were receding from the 
Lya source too rapidly (see § 15. 1[) . 
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Figure 7. The left panel shows the predicted absorption line strength as a function of impact parameter for model IV, while the right 
panel shows the predicted Lya surface brightness profiles as obtained from the Monte-Carlo simulation {black filled circles), and analytic 
calculation (black dashed line). This model - in which gravity decelerates the cold clumps after their initial acceleration - clearly predicts 
a surface brightness profile that is much closer to the observations than model I-III. 
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Figure 6. This plot shows the velocity profiles that we adopt for 
model I {solid line), and model IV {dot-dashed line). In contrast 
to model I, in model IV gravity causes the clumps to decelerate 
beyond some critical radius r^rit- As a result, clumps flow out at 
lower velocities in model IV, which makes them more opaque to 
Lya photons, which results in brighter Lya halos (text). 



Tile goal of this section is to investigate wlietlier we 
can reproduce tlie absorption line and Lya lialo data better 
in simplified models wiiicii take tliis gravitational decelera- 
tion of tiie cold clumps into account. Following Murray et 
al. (2005, also see Martin 2005) we write tlie momentum 
equation for a cold, optically thick, clump as 



dvc 
~dt 



GM{r) 



+ Ar 



(13) 



where we used the same function to describe the cloud ac- 
celeration as before fS 13 1[1. The case a = 2 , A = 2a^Rg 
corresponds to Eq 24 of iMurrav et all l|2005l ). Here, Rg = 
rinin{L/Lcdd), where the ratio L/L^dd denotes the total lu- 
minosity of the source normalized to the Eddington luminos- 
ity of the galaxy ( s ee Mu rray et al. 2005 for details). Follow- 
ing [Murray_eLalJ (120051) we take for the model of the grav- 



itational potential that of an isothermal sphere for which 
M(r) = 2a^r/G, in which a denotes the velocity dispersion. 



The solution to Eq[T3]is given by 



Vc{r) = 2a 



f(¥ 



-f 



A 



2cr2(l 



(14) 



where we assumed the boundary condition Vc{rinin) = 0. 

In theory, it is straightforward to repeat the analysis of 
§ 13.11 and do a MCMC simulation to simultaneously con- 
strain all model parameters including rmin. A, and a by 
finding the best-fit model to the absorption line data. We 
have instead fixed the values a = 1.4, rmin = 1 kpc, which 
we found to provide good fits for model I-model III. We 
also assumed a = 150 km s~^, which is the value that we 
assumed previously (see §[5]), and -Rg = 2.5 (i.e. A — 5(t^). 
This latter choice is entirely empirical: The dot-dashed line 
in Figure [6] shows that the resulting velocity profile fEQ ll4|) 
reaches a maximum of «c,max ~ 350 km s~^ at r ~ 10 kpc 
and then decreases to Vc ~ 100 km s~^ at r = 250 kpc. For 
smaller values of -Rg, the clumps would turn around and fall 
back onto the galaxy, while larger values of -Rg would result 
in negligible deceleration. The solid line shows the velocity 
profile that was adopted in model I. 

The MCMC finds the best-fit model for the parameters 
shown in Table [51 This model - which we refer to as model 
IV - is compared to the absorption line data in the left panel 
of Figure [71 The right panel of Figure [71 shows the predicted 
Lya surface brightness and polarization profiles. It is clear 
that the predicted surface brightness is significantly higher 
than that predicted by model I-model III and agrees with 
the data to within a factor of ~ 2 — 3. The enhanced surface 
brightness profile is a direct result of the lower outfiow ve- 
locity of the clumps which makes them more opaque to Lya 
photons. The photons still most often scatter off 1 clump, 
and the analytic calculation still compares quite well to the 
result we obtained from the Monte-Carlo simulation. The 
predicted polarization is also high for this model: with the 
linear polarization reaching P ~ 40% at a surface brightness 
level of SB~ 10~^* erg s~^ cm~'^ arcsec"^. 
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Figure 8. This Figure shows surface brightness contours for a biconical outflow model as seen along the axis of the cone {left panel), 
an perpendicular to the cone axis {right panel). The black/ red/ blue contour denote a surface brightness level of log[SB/(erg s~^ cin~-^ 
arcsec~^)]=-18.0/-19.0/-20.0. The opening angle of the cone is 90°. This plot illustrates that a case in which our code works on non- 
spherical clump distributions. 



Importantly, this model still predicts that /ns ~ 60% of 
the photons does not scatter at all in the outflow. This model 
therefore also predicts that the Lya halo is accompanied 
by a point source of comparable luminosity, which is not 
observed. 



6.2 Model V: Concentrating the Outflow into 
Cones 

The flnal model that we consider - which we refer to as 
model V - is a biconical outflow model, in which the hot 
and cold gas escape from the galaxy along two cones whose 
axes are parallel. The reason that we study such a model is 
the observed azimuthal dependence of the Mg II at 6 < 50 
kpc of inclined disk-dominated galaxies at z = 0.5 — 0.9, 
which indicates the presence of a biconica l outflows that ar e 
aligned along the disk rotation axis (Bordoloi et al.ll2011f ). 
We would like to see how our results change if we drop the 
assumption of having symmetric outflows. 

We take the parameters from model IV but now com- 
press the outflowing clumps into two cones. For simplicity 
we align the cone axes with the z-axis, and assume that the 
opening angle is 9 = 45° (i.e. the edge of the cone intersect 
the z-axis at 45° ) . This compression of the outflow has two 
implications: (j) the number density of cold clumps is en- 
hanced by a factor of 1/[1 — cos 9] ~ 3.5 at each radius, and 
{a) the pressure of the hot wind also increases by a factor 
of 3.5, which compresses the radius of the cold clouds by a 
factor of 3.5^'^ = 1.5. The dotted Ime in the left panel of 
Figure [7] shows the predicted EW vs b curve for a random 
realization of model V. This model also nicely reproduces 
the observations without any further tuning. We have not 
performed an analytic calculation of EW as a function of 6 
(as in EqO, as this would require averaging over random 
cone orientations. 

Figure [S] shows surface brightness contours for model 
V when we view the outflow along the cone axis {left 
panel), and perpendicular to the cone axis {right panel). 
The black/ red/ blue contour indicates a surface brightness 
level of log[SB/(erg s~^ cm'^ arcsec-2)]= -18.0/-19.0/-20.0. 
This Figure illustrates that our code works on non-spherical 



clump distributions. For this calculation, we assume that all 
Lya photons escapes from the galaxy into the cones (i.e. the 
central source does not emit Lya photons isotropi- 
cally, but instead into two cones w^ith opening angles 

if 45°), which is motivated by the physical picture in which 
the outflow has 'cleared' out a cone of lower H I column 
density along which Lya photons escape. Studies of Lya 
transfer through simulated galaxies have also shown that 
the escape of Lya photons from the dusty ISM of galaxies 
can proceed highly anisotropically: the Lya flux transmitted 
i nto different directi o ns may va ry b y an o rder of magnitude 
(Laursen et al., ,200a : ,Yaiima et al. '.2012bh. 

The right panel of Figure [7] show the predicted 
azimuthally-averaged surface brightness & polarization pro- 
files when the bipolar outfiow is viewed from the top {open 
circles), and from the side {open triangles). The predicted 
polarization is lower (higher) when we view the outfiow from 
the top (side), as photons typically scatter hy fi < cos 45° 
{fi > cos 45°). This panel also shows that the azimuthally 
averaged surface brightness profile depends only weakly on 
whether the outfiow is bipolar or not, and on from which 
angle we view the outfiow. This last result is important as it 
can alleviate the problem that model IV has, namely that 
~ 50% of all Lya photons do not scatter in the outflow, 
and should thus be detectable as a Lya source. Invoking 
this bipolarity may help us avoid predicting a Lya point 
source with Lya halos, because when 9 = 45°, a fraction 
cos6' ~ 0.71 of all the sightlines would not lie in the cone, 
and we would not see a point source along these sightlines. 

Two caveats to model V are: {i) as mentioned and jus- 
tified above, our model assumes that Lya photons escape 
from the galaxy into the cones of outflowing gas. This 'fo- 
cussing' of Lya photons into the biconical outflow requires 
additional H I gas which is not present in our model. This 
gas would contribute to the measured EW at 6 = 0, and 
would hence reduce the amount of H I gas that we are al- 
lowed to assign to the clumps at small impact parameters. 
This reduced amount of cold gas would suppress the pre- 
dicted amount of scattering - and hence the predicted sur- 
face brightness - at small impact parameters (which would 
agree better with the data); {ii) if the outflowing material 
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were indeed all confided to two cones with opening angles 
9 = 45°, then we fail to explain why blueshifted absorption 
lines of low ionization abs orption lines are observed in prac- 
ticall y all sightlines fe.g. IShaplev et aljboosl : ISteidel et al.1 
I2OI0I ). We intend to include the additional constraints pro- 
vided by the 'down-the-barrel spectra' in a more realistic 
model of the outflows, in which the clump distribution is not 
only a function of radius r, but also of azimuthal angle (see 
§[Z1I- 



7 DISCUSSION 

T.l Discussion of Model Uncertainties &: Caveats 

Here, we discuss potential implications of our adopted as- 
sumptions and simplifications. 

• We constrain the HI content of the cold clumps in the 
large-scale outflow with the mean observed absorption line 
strength (EW) as a function of impact parameter (fe). These 
same clumps scatter Lya photons emitted by the foreground 
galaxy into the line-of-sight, which could weaken the ob- 
served line strength (see Prochaska et al. 2011). For example, 
the observed surface brightness at & = 31 kpc is 6 x 10~^^ 
erg s~^ cm~^ arcsec"'^. The Lya absorption strength has 
been measured over an area of 1.2 x 1.35 arcsec^ (C. Steidel, 
private communication), which implies that the observed ab- 



10" 



erg s cm of scat- 



sorption is accompanied by 
tered radiation. 

We can compare this scattered flux to the fiux that has 
been absorbed out of the spectrum of the background galaxy 
as follows: (i) the observed fiux density (/a in erg s~^ 
A~^) at restframe A = 1100 A is comparable to that at 
A = 1500 A (see Fig 5 of Steidel et al. 2010, and using that 
/^ ex A^/a); (ii) we estimate the observed restframe fiux 
density at A = 1500 A (restframe) from the average UV- 
derived star formation rate, uncorrected for dust, which was 
SFR= 8 Mq yr~^ assuming the K ennicutt (1998) star for- 
mation calibrator (jErb et al.ll2006l ). For a galaxy at z = 2.3 
this star formation rate translates to an observed restframe 
/a(A = 1500 A) ~ 2 X 10"^* erg s"^ cm"^ A'^ After 
combining (i) and (ii), we find that the scattered fiux of 
~ ]^Q-is gj.g g-i c]-n-2 corresponds to a rest frame EW~ 0.5 
A, which lies a factor of 4 below the measured EW. This 
implies that scattered flux is sub-dominant to the absorbed 
flux. At larger impact parameters the scattered flux becomes 
exponentially fainter and even less important. 

Obviously, our calculation is approximate and relies on 
an average spectrum, an average star formation rate, and 
a median redshift. In certain cases, we expect the scattered 
flux to be more important than our previous estimate. How- 
ever, in such cases the spatially extended scattered flux is 
detectable in spectra of pixels adjacent to the background 
galaxy, and can be corrected for. We therefore conclude that 
the potential contribution of scattered radiation to the ob- 
served EW does not affect our results. 

• When constraining the HI content of the clumps (using 
the observed EW as a function of 6), we implicitly assume 
that only the cold clum ps in the outflow contribu te to the 
observed EW. However, Ivan de Voort et al.l (|201l|) have re- 
cently found that in their hydrodynamical simulations, the 
cold streams of gas that are feeding the central galaxy may 



produce more large column density absorbers (A^hi JjlO^" 
cm~^) around galaxies than outflows (also see Funmgali et 
al. 2011, but note that modeling the outflowing component 
at large impact parameter is extremely uncertain, see §[3). 
The probable contribution of cold streams to the observed 
EW at a given impact parameter reduces the amount of HI 
that we can assign to cold clumps in the outflow. With less 
HI in the outfiowing clumps, we expect them to scatter fewer 
Lya photons, and that consequently our surface brightness 
levels are reduced by an amount which depends on the the 
overall contribution of the outflowing cold clumps to the 
observed EW at a given b. 

• The observed EW as a function of b has b een measured 
aroun d galajcies with a mean redshift (2) — 2.2 (jSteidel et al.l 
120101 ). while the mean redshift of galaxies that were used in 
the stack of narrowband images was (zhaio) = 2.65. While 
galaxies in both samples were both selected in very simi- 
lar ways, and therefore likely trace similar populations (e.g. 
both populations have virtually identical mean star forma- 
tion rates, see § 13. 1|) . the observed Lya absorption line 
strength does not really probe the same gas that is scatter- 
ing the Lya halos. For this reason, we consider differences 
between the predicted and observed Lya surface brightness 
profiles at the factor of 2-3 level (as was the case in model 
IV) not a problem. On the other hand, it is unrealistic to at- 
tribute the order of magnitude differences in the predicted 
surface brightness profiles (which we found for model I- 
model III) to selection effects. 

• In our model in which the clumps decelerate after their 
initial acceleration (see §[67l]), the clumps reach a maximum 
velocity of «c,max ~ 350 km s~^, which is well below the 
maximum velocity that was inferred by Steidel et al. (2010) 
of Umax ~ 800 km s~^. Fujita et al. (2009) found in their hy- 
drodynamical simulations that a small fraction of the cold 
shell fragments could be accelerated to reach large outflow 
velocities ( ^750 km s~^), while the bulk of the gas was trav- 
eling at lower velocities of ~ 300 km s~^. This suggests that 
the maximum velocity that is inferred from the observations 
can be consistent with our model which only describes the 
kinematics of this bulk of the gas. Furthermore, because of 
the large HI column densities in the clumps, the photons 
can scatter in the wings of the Lya line proflle. As a result, 
Lya absorption may trace a wider range of velocities than 
the range of actual outflow velocities. 

• In our model there is a one-to-one mapping between 
radius and velocity. In reality, we expect outflows to have 
a range of ve locities at a given radius: e.g. in the model of 
iMartinI (|2005|), the cloud acceleration increases with their HI 
column density. For a range of HI column densities we there- 
fore expect the cold clouds to have a range of velocities at a 
given radius. This scatter may boost the predicted Lya sur- 
face brightness profile, in particular when this scatter gives 
rise to a population of clumps with a lower Vc{r). 

• Our model assumes that there is a unique clump mass. 
In reality, there is a distribution. This is very likely not an 



10* M( 



0) 



issue: in model I the clump mass was rriciuui 
while in model II the clump mass was ~ 10 times larger. 
As long as the distribution of clumps is constrained by the 
absorption line data, we predicted virtually identical surface 
brightness profiles. We therefore consider it unlikely that we 
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obtain significantly different surface briglitness profile, if we 
assume a finite range of clump masses. 

• Our model assumes a Gaussian emission profile, for 
which the width is set by the velocity dispersion of the gas. 
However, the outflow contains swept-up shells of cold neu- 
tral gas before it breaks out of the galaxy. It is therefore not 
unlikely that Lya photons scatter off these dense cold shells, 
which would result in an overall redshift of the line, before 
the Lya photons escape from the galaxy mto the large scale 
outflow. Redshifting of the Lya line would reduce the overall 
scattering probability in the outflow, and could reduce the 
predicted surface brightness profiles. 

• We assumed that the escape fractiorO of Lya pho- 
tons was 50-100% (depending on the choice of IMF, and 
gas metallicity as these affect the intrinsic Lya luminosity 
of the galaxy at a fixed SFR) . This potentia l ly hig h escape 
fraction was already noted by ISteidel et al.l (|201lh . In our 
models, it is required to be higher by a factor of ~ 2 be- 
cause about half of the photons that escape from the galaxy 
into the large scale outflow were not scattered at all, and are 
effectively wasted. 

This radiation that is not scattered in the outflow must be 
a point source of comparable (model IV) or higher (model 
I -model III) luminosity than the halo itself, which is in 
conflict with observations. This disagreement can be alle- 
viated by invoking that the outflow is bipolar (see § 16. 2|) . 
and/or by a population of low column density absorbers (as 
observed in galaxy-quasar pair data by Rudie et al. 2012) 
that have a velocity Vc that differs substantially from that 
in of the clumps in model I-V (see § 17. 3p . 



7.2 Connection with Lya Blobs 

ISteidel et all |20lJ) found that the Lya 'blobs' - defined 
loosely as having an area ^16 arcsec"^ for which SB J>10~'^* 
erg s~^ cm~^ arscec"^ (Matsuda et al. 2004, also see Steidel 
et al. 2000, Saito et al. 2006) - within their narrowband sur- 
vey volume to have surface brightness profiles almost identi- 
cal to that of the Lya halos, apart from an overall off-set in 
the surface brightness. Scattering of Lya photons in a large 
scale outflow can theoretically explain Lya blobs, if we in- 
creEise the star formation rate of the central Lya source. This 
is illustrated in Figure |9] which shows the dependence of the 
predicted surface brightness profile on SFR for model IV 
(which was also shown in Fig [7]) . In these calculations we 
kept all other model parameters fixed. 

The dashed line shows the original prediction for model 
IV which assumed SFR=34 Mq yr"^ The dotted line/ dot- 
dashed line shows the predicted surface brightness profile for 



SFR=100 Mq yr-i / SFR=10 M, 



yr 



The predicted sur- 



face brightness profile depends quite strongly on SFR. This 
is because: («) the intrinsic Lya luminosity of the central 



^■^ This escape fraction refers to the fraction of photons that es- 
cape from the ISM of the galaxy into the large-scale outflow. 
This escape fraction differs from the e scape fraction that is used 
in recent observa tional papers (e.g. iHaves et al.l |2010| . l2011bl : 
iBlanc et al.ll201ll) . which represents the ratio of the observed to 
the intrinsic Lya flux, which can depend sen sitively on the sur- 
face brightness threshold of the observations llZheng et al.ll201CI . 
also see Jeeson-Daniel et al. in prep). 
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Figure 9. This plot shows the dependence of our surface bright- 
ness profiles of the Lya halos, as a function of the star formation 
rate (SFR, in Mq yr~^), assuming the same model parameters 
for the outflow (mass loading factor, number of clumps, etc.). The 
predicted surface brightness profile depends strongly on SFR be- 
cause (i) the intrinsic Lyo luminosity of the central source scales 
linearly with SFR, and (m) the total amount of scattering ma- 
terial scales linearly with SFR. As a result, changing SFR by a 
factor of ~ 3 can change the surface brightness at a fixed B by an 
order of magnitude 



source scales linearly with SFR, and [ii) the amount of cold 
gas that can scatter the Lya photons also depends linearly 
on SFR. These two effects combined suggest that the surface 
brightness at a given impact parameter can depend on SFR 
as oc SFR^, which can explain that the surface brightness 
at a given impact parameter can vary by an order of magni- 
tude as a result of a factor of ~ 3 change on the SFR (which 
gives more weight to the contribution of high-SFR galaxies 
to the stacked Lya image). The feature in the SFR=10 Mq 
yr^^ curve at S ~ 9 arcsec refiects that the clumps do not 
self-shield at r ^76 kpc in this model. 

Figure [9] shows that it is possible to have Lya 
blobs around galaxies that are forming stars at a rate 
SFR ^lOOAf© yr~^, as the predicted surface brightness 
profile drops below ~ 10~^* erg s~^ cm~^ arcsec"^ only 
at 6 ^7 — 8 arcsec. There is observational evidence that 
a fraction of the Lya blobs are associated with sub-mm 
galaxies, which are bel ieved to be forming stars at rates 
of SFR^ I Q ^Mq yr~^ JChapman et al.ll200ll: iGeach et all 
l2005l . I2OO7I : iMatsuda et all l2007f ). IChapman et all (|2005l ) 
have noted that a remarkably high fraction (~ 50%) of sub- 



m m galaxies show Lya em ission lines in their spectra (also 
see iNilsson fc M0lleil |2009| . for a detection of Lya from a 
ULIRG). This may indeed suggest that enough Lya pho- 
tons escape from the dusty interstellar media of sub-mm 
galaxies, and then scatter in the large scale outfiow to ac- 
count for the blobs. Support for the notion that Lya blobs 
consist of scattered radi ation is pro vided by the detection of 
polarization by Hayes e t al.l (|2011al ). albeit at a level that is 
lower by a fact or of ~ 2 than the val ues predicted here. On 
the other hand. lPrescott et al. 11201 1') put an upper limit on 
P ^10% in LABdOS (JDev et al.l 12005'). which clearly rules 
out our models. 

It appears increasingly plausible that there are distinct 
physical mechanisms that power Lya blobs. The polariza- 



© 2006 RAS, MNRAS 000.[lHT9l 



16 Dijkstra & Kramer 



tion measurement of lHaves et alj (l2011al ) clea rly favors mod- 
els tha t invoke scattering. On the other hand. lPrescott et alJ 
((2013) find that the UV continuum (non-ionizing) associated 
with LABd05 is also spatially extended, which favors having 
spatially extended emission of both UV and Lya photons. 
This observation could be consistent with dust scattering, 
which would explain why the polarization of Lya would not 
have been detected (see §|47l]). However, given the large ob- 
served EW of the diffuse Lya when m easured relative to th e 
diffuse UV-continuum (REW» 200 A. IPrescott et al.ll2012l ). 
this is not very plausible. Alternative support for the notion 
that some blobs are not a result of scattering is provided by 
thos e blobs that do not have any clear galaxy counterparts 
(e.g. iNilsson et ai]|2006l : ISmith fc Jarld?l200'7l ) . 



7.3 Additional Constrains from Galaxy-Quasar 
Pair Data 



7.4 Outlook & Potential Improvements 



iRakic et al] l|201lf ) and iRudie et all (|2012l l use galaxy- 
qua sar pairs to probe the CGM of the foreground galax- 
ies. iRudie et al.l (|2012f ) present 10 pairs for which the back- 
ground quasar lies at b < 100 kpc. They find 6 absorbers 
with log A''hi ^17.0, of which 2 absorbers have log A'^ei ^18.0, 
of which one absorber has log A^hi ~ 20 . We compare this to 
model V at fe = 71 kpc, which corresponds to the median 
value for b if the sightlines are distributed randomly within 
the circle for radius b — 100 kpc. We find that our model pre- 
dicts ~ 7 absorbers with log Ahi ^17.0, of which ~ 7 have 
log Am ^18.0, of which ^1 absorber has log Am ^19.3. 
Given the simplified nature of our model, and the relatively 
small number of observed sightlines, we consider these num- 
bers encouraging. For example, our predicted number of ab- 
sorbers with log Ahi >18.0 is reduced to ~ 3, and therefore 
more consistent with observations, if we simply increase the 
critical number density above which gas self-shields by a 
factor of 2 to ncrit = 0.012 cm~^, which is still reasonable 
(e.g. Faucher-Giguere et al. 2010 adopt ncrit = 0.01 cm~^). 
Such a modification reduces the predicted EW as a function 
of b, but only significantly at b ^70 kpc. In this case the 
observed EW at 6 = 100 kpc could be accounted for by a 
large number of lower column density absorbers, which have 
been observed (see below) but which are not present in the 
model. 

Another difference worth emphasizing is that our model 
predicts that ten sightlines wit h b = 71 kpc should intersect 
a total of ~ 13 cold clumps. iRudie et al.l (|2012l ) find sig- 
nificantly more low column density, log Ahi = 14.5 — 17.0, 
absorbers. This implies that our clumpy outflow model does 
not account for all the observed absorbers, and thus all po- 
tential 'scatterers'. However, if we wish to use low column 
density absorbers to scatter photons into Lya halos, then 
they must have lower outflow velocities (or they have to be 
inflowing) than the clumps in our model, otherwise they are 
transparent to the Lya photons. The possible presence of 
these low-column density that move at different velocities 
than the high column density clumps in our models, may 
have the interesting benefit that they reduce the fraction 
of photons that do not scatter in the outflow at all. These 
clumps may thus reduce the luminosity of the central Lya 
sources that accompanies the Lya halos in model I-IV, 
and in model V when we view the outflow down one of the 



In § 17.11 we highlighted the simplifications of our model, 
which underlined that many improvements are possible. We 
discuss some examples of how we intend to improve upon 
our analysis below. 

We have so far focused on using the observed Lya ab- 
sorption line strengths as well as the surface brightness pro- 
file of Lya halos to constrain parameters of outflow models. 
However, there is information encoded in the observed spec- 
tral l ine shape of the Lya emission line (e.g. lYamada et al.l 
|2012| . and references therein), as well as whether it is red- 
shifted or blueshifted relative to the galaxies systemic veloc- 
ity. For example, Lya lines that are blueshifted/redshifted 
with respect to the systemic velocity - to first order 
- are indicative of scatteri ng through an opaque inflow- 
ing/outflowing medi um (e.g. lZheng fc Miralda-Escudell2002l : 
IPiikstra et al.l l2006h. The first joint Ha - Lya s pectral line 
obser vations of spatially extended Lya nebulae (|Yang et al.l 
l201ll l. compact Lya selected galaxies (Finkelstein et al. 
2011, also see McLinden et a 2011 for joint Lya-[0 III] ob- 
servation s), and 'double-pea ked' Lya emitting UV-selected 
galaxies (IKulas et al.l 120121 ) have recently been reported. 
iKulas et al.l ([2013) have already shown that such observa- 
tions can rule out the 'shell-models' for outflows. In the fu- 
ture we plan to explore what additional constraints we can 
place on outflow models on small scales (i.e. r _^10 kpc) 
with observations of the Lya spectral line shape (and shift). 
Finally, our previous discussi on (§17.311 showed t hat observa- 
tions of galaxy-quasar pairs (JRakic et al.ll201ll : iRudie et al.l 
120121 ) already provide useful additional constraints on our 
models. 

We also plan to extend our stu dy to include other lines. 
For example, iBordoloi et al.l ()201lh have presented the radi- 
ally (and azimuthally) dependent absorption line strength of 
Mg II around bright flux-selected galaxies at 0.5 < 2 < 0.9 
from the z- Cosmos redshift survey. We can constrain the 
Mg II content of the clumps in our model by matching this 
data. We can then make predictions for surface brightness 
profiles of scattered Mg II emission, and compare to ob- 
servations of spatially ext ended Mg II emis sion around a 
z = 0.69 starburst galaxy l|Rubin et al.ll20lj) . 



8 CONCLUSIONS 

We have presented 'constrained' radiative transfer calcula- 
tions of Lya photons propagating through clumpy, dusty, 
large scale outflows, and explore whether scattering through 
such an outflow can quantitatively explain the Lya ha- 
los that have been observed around Lyman Break Galax- 
ies (LBGs, see Steidel et al 2011). As part of our analy- 
sis we have modified a Lya Monte-Carlo radiative transfer 
code to allow us to follow the propagation of Lya photons 
through a multiphase, dusty medium for arbitrary distribu- 
tions of clumps. This code also computes the polarization 
of the scattered Lya radiation. Previous calculations of the 
polarization of scattered Lya radiation have focused only 
on homogeneous spherically symmetric gas clouds or shells. 
We have successfully tested our code against several ana- 
lytic solutions, some of which - in particular the directional 
dependent frequency redistribution function - are new. 
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Modeling the distribution and kinematics of cold gas in 
outflows from first principles is an extremely complex task, 
which likely requires magneto-hydrodynamical simulations 
that have sub-pc resolution (§ [3}. We have taken an dif- 
ferent approach, and constructed phenomenological models 
for the large-scale outflows in which cold (logTc ~ 3 — 4) 
clumps are in pressure equilibrium with a hot (logTh ~ 7) 
wind. We first considered models in which the cold clumps 
are distributed symmetrically around the source, and which 
accelerate continuousl y as they break out of the interstellar 
medium of the galaxy. Steidel et al.l ([2010) showed that this 
type of model may qualitatively simultaneously explain the 
observed Lya absorption line strength in the CGM, as well 
as the observed surface brightness profiles of Lya emission 
line halos. 

Our more detailed analysis shows that such models- 
which contain 10^~^ discrete clumps-can reproduce the 
observed Lya absorption strength of the circumgalactic 
medium measured in the spectra of background galaxies very 
well, and for model parameters that are physically plausi- 
ble (see Fig[T|. However, when we insert a Lya source in 
the center of these clumpy outflow models, and compute 
the observable properties of the scattered Lya radiation, we 
typically flnd that the predicted Lya halos are significantly 
fainter and more concentrated than what is observed (Fig|4]). 
The reason for this discrepancy is easy to understand: out- 
fiowing cold clumps that scatter photons at large (6 ^30 
kpc) impact parameters are propagating away from the Lya 
source at v ^600 km s"'^. In order for the clumps to scatter 
Lya photons they must be opaque to these photons, which 
requires HI column densities in excess of A^hi JjlO^^ cm"'^. 
However, the absorption line data requires that the num- 
ber of such clumps at large impact parameters is so small 
that a significant fraction of the photons never encounter 
them. Our conclusion that we cannot simultaneously fit the 
absorption line and Lya halo data -with a clumpy outflow 
in which the clumps are distributed spherically around the 
galaxy and which accelerate with radius - is therefore ro- 
bust. 

We also found that the vast majority of photons scatter 
in zero or one clumps, and that it is possible to analytically 
compute the Lya surface brightness and polarization pro- 
files (see Fig |4)| . The fact that a significant fraction of the 
photons do not scatter in the outflow is problematic. The 
photons that do not scatter must be detectable as a point 
source, and it's predicted luminosity equals or exceeds the 
total luminosity of the Lya halo, which is in further dis- 
agreement with the observations. 

We can much better simultaneously reproduce the ob- 
served Lya absorption line strengths and the Lya halos 
with models in which the cold outflowing clumps decelerate 
(see Fig [7]). This deceleration occurs naturally in models of 
momentum-driven winds (Jj IG.lfl . We can alleviate the prob- 
lem of predicting a bright Lya point source to accompany 
the halo, if the outflow is bipolar with an (half) opening 
angle 9 ;^45° (i; l6.2[) . This problem may be further reduced 
if the observed additional low-column density absorbers in 
galaxy-quasar pairs move at velocities that allow them to 
resonantly scatter an additional fraction of the Lya photons 
(see §[73J. 

We found that models which do flt both the absorption 
line strength and Lya halo data give rise to levels of linear 



polarization that reach P ~ 40% at a surface brightness level 
of SB~ 10"^® erg s"^ cm"^ arscec"^ (e.g. Fig[7l). This po- 
larization signature is likely unique to the scattering models 
and likely distinguishes them from models in which the Lya 
photons were emitted over a spatially extended region (see 
§ [T} . Furthermore, because the large polarization signature 
is a result of non-resonant scattering, the polariza t ion als o 
distinguishes our models from those of lZheng et al.l (l2011al ). 
in which the halos were a result of resonant scattering (§[T|. 
It should be noted that it remains to be shown that pre- 
dictions of these other models are in quantitative agreement 
with the observed Lya absorption line data, as well as the 
Lya halos. 

This paper illustrates clearly that Lya emission line ha- 
los around star forming galaxies provide valuable constraints 
on the cold gas distribution & kinematics in their circum- 
galactic medium, and that these constraints nicely comple- 
ment those obtained from absorption line studies alone. 
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APPENDIX A: CODE DESCRIPTION & 
TESTING 



We d escribe modifications to the code of iDiikstra et al.l 
1(2003) ill detail. We test various subroutines of the new 
code in a simple geometry in which six clumps of radius 
Re ~ 10 kpc lie on the coordinate axes a distance d = 50 
kpc from the origin. That is, the clumps lie at {x = ±d, 0, 0), 
(0, y — ±d, 0), and (0, 0,z — ±d), where d = 50 kpc. We as- 
sign an outflow velocity Vc to the clumps 

In section § lAll we consider cases in which the clumps 
are transparent (i.e. t ^ 1) to Lya, and in §[X2]we consider 
cases in which the clumps are extremely opaque to Lya pho- 
tons. For all these tests we assume that the gas temperature 
of the gas in the clumps is Tc — 10^ K. 




y 



Figure Al. Clump geometry that we assumed in our test cal- 
culations. Six clumps of radius Re lie on the coordinate axes a 
distance d = 50 kpc from the origin, which contains the source of 
Lya photons. When the clumps are optically thin to Lya photons, 
we can analytically compute the fraction of photons that scatter 
in N clumps as a function of A^ (see § lAll l . When the clumps are 
opaque to Lya photons (see § IA2I I , we can analytically compute 
the emerging spectrum. 



Al.l Testing the Interclump Propagation Scheme 

The sky covering factor of a single clump for a central source 
is /c = 4^-, where Q,c — irO^ in which 9c = arcsin(_Rc/d). 
The fraction of photons that scatter in only one clumjo is 
/iciump = 6x/cX (l-exp[-r(a;in)]). Here {l-exp[-r(a;in)]) 
is given by 



27r 



ydy 1 - exp 



To X 



i?2 



exphr(xin)]) = (Al) 



R^ 



Here y denotes the impact parameter from the center of the 
clump to where the photon strikes, to denotes the line center 
optical depth through the center of the clump, and <j){^'[y]) 

denotes the Voigt profile evaluated at x' = Xi^ ^^ cos 61, 

where 9 denotes the angle between the photon's wavevector 
k and the outflow velocity vector v, and we have cos 9 = 
\/l-sin2 0= Vl-(y/d)2. 

For example, for stationary clumps, the y— dependence 
of the term (f>{x'[y]) vanishes and the integral can be eval- 
uated analytically when to <^ 1, resulting in /idump ~ 
0.0405To(?i(a;). The left panel of Figure IX2I compares analyti- 
cally computed values of /idump for Wc = {solid line), and 
Vc = 20 km s~^ {dotted lines) as a function of to, with those 
obtained from the Monte-Garlo code. The agreement is ex- 
cellent, and shows that inter and intra-clump propagation 
are captured well by the Monte-Garlo code. 

In the optically thin limit, we can also estimate the 



Al Central Source & Transparent Clumps 

For these tests, we insert all Lya photons at the origin, and 
at the line resonance, i.e. a;in = 0.0. 



^^ Formally, we have to multiply this probability /iscat by the 
probability that photons do not subsequently scatter in other 
clumps. As the clump sky covering factor is only ~ 6%, this in- 
troduces a correction of at most a factor of ~ (1.0 — 0.06) = 0.94. 
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Figure A2. Comparison between analytic and Monte-Carlo calculations of the fraction of photons that scatter in one clump (left panel) 
and in two distinct clumps [right panel) as a function of the line center optical depth of an individual clump, tq, for Vc = (solid lines), 
and Vc = 20 km s~^ (dotted lines). The agreement is excellent, and shows that inter and intra-clump propagation-as well as frequency 
redistribution by moving clumps— are captured well by the Monte-Carlo code. 



fraction of photons that scatter in two different clumps. An 
accurate estimate for this probability is given by 



!^c,i 



/2clump — /iclump ^ / ^ . -* vM 



«#j 



I dx' \{1 - exp-Ti{x'))]R(x'\x[xin], m). (A2) 

J — oo 

This equation gives the probability that a photon scatters 
in a second clump, denoted with number 'i', after having 
scattered in the first clump, denoted with 'j'. The probability 
that the photon scatters in a second clump is a product of 
the probabilities that the photon scatters into a sightline 
that intersects clump 'i', and that it then scatters in that 
clump. This latter probability depends on the frequency of 
the photon after the first scattering event, and we integrate 
over all possible photon frequencies weighted by PDF of this 
frequency. A more detailed quantitative explanation follows 
below. 

Firstly, the probability that the photon scatters into a 
sightline that intersects clump 'i' is given by ~ -^P(fii). 
In our testcase, a photon has to scatter either by /i ~ —1 
for scattering in the clump on the same coordinate axis (i = 
5), or /K « —1/^/2 for scattering in clumps on one of the 
other coordinate axes (i = 1 — 4). These probabilities are 
approximations -but accurate ones- because in reality the 
photons scatter into a (narrow) range of fi, which in detail 
depends on where exactly the first scattering event occurred. 
To capture this effect properly, we would have to average 
over /i weighted by the proper PDF for /i. However, since this 
range of fi only extends over A/i ~ 0.2, and because both 
P(l-i) and R(x'\x[xin], fJ.i) change very little over this range, 
this more detailed and tedious procedure barely changes our 
final results. 

Secondly, the expression for the probability that a pho- 
ton scatters in the second clump is given by (1— exp —Ti(x')), 
which is given by Eg I A2I for i — 5 (with d = 100 kpc), but 
for i = 1 — 4 we omit the y— dependence of x' . The geometry 
for scattering in clumps 1 = 1 — 4 does not allow for a simple 



mapping between impact parameter y and Doppler boost, 
and this last modification represents a reasonable approxi- 
mation. 

Finally, the PDF for the outgoing photon frequency x' 
depends on both the scattering direction and incoming pho- 
ton frequency iin. We derive an analytic solution for the 
frequency redistribution function, R(x"\x, n), which denote 
the PDF for x" given x. These frequencies are measured 
in the frame of the gas. In our test case, photons appear 
at frequency x = Kin — ^c/uth in the frame of the first 
clump. The outgoing photon frequency x' (measured in the 
lab frame) relates to x" through a Lorentz transformation: 
x' — x" + fiVc/vth- The expression for R(x"\x, n) is derived 
in Appendix [Cl 

The right panel of Figure IX2l compares the analytically 
computed values of /2ciump as a function of tq for Uc = 
(solid line), and Vc = 20 km s~^ (dotted lines) with those ob- 
tained from our Monte-Carlo code. The agreement is again 
excellent. This further illustrates that inter and intra-clump 
propagation of photons is described accurately. Further- 
more, the (directional dependent) frequency redistribution 
is also captured accurately. 



A 1.2 Testing the Surface Brightness & Polarization 
Subroutines 



iDiikstra et al.l (|200q ) focused on spherically symmetric gas 
distributions, which makes the calculation of the predicted 
surface brightness distribution straightforward. In three- 
dimensional geometries ideally one has to compute the frac- 
tion of photons that escapes from the medium, exactly in the 
direction of the telescope. Since the telescope is a cosmologi- 
cal distance removed from the location of last scattering, and 
hence practically subtends an infinitesimally small fraction 
of the sky, this procedure is not possible in practice. 

We follow a standard approach for computing surface 
brightness profiles, and compute the differential probabil- 
ity that a photon is scattered exactly towards the tele- 
scope, for each scattering event. This probability can be 
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Figure A3. The surface brightness profiles for the scattered radiation are shown for all 6 viewing directions for a case in which tq = 0.1. 
The lower right panel shows the average surface brightness profile. The central clump is approximately twice as bright as this contains 
emission that scattered from two clumps. This plot illustrates nicely that all scattered radiation is confined to regions associated with 
clumps. A more detailed view of a the surface brightness profile of a single clump in the averaged image is shown in Figure [A4l 



computed as follows. A photon scatters at some location 
Xscat = {xs,ys, Za), by an atom whose velocity components 
are given by v = (wx, Uy, Wz), or by a dust grain whose ther- 
mal motion we neglect. We denote the photon's propagation, 
frequency, and polarization before scattering with kin, xin, 
and Gin. Now let us consider the a:"*" — y^ image. Photons 
that make up this image would have to be scattered into 
direction kout = (fcx, fcy, fe) = (0, 0, 1). 

The probability per sterradian that a photon escapes 
into direction kout is 



P = exp[— T(3;out, kout, Xout)] x 
xP(kin, kout, ein I wing/res/dust). 



(A3) 



where a;out is determined fully by fcout, a;in and the veloc- 
ity vector of the scattering atom, v (see Ea lCl|l . and where 
P(kin, kout,ein|wing/res/dust) denotes the phase function, 
which depends on whether the photon scatters in the line 
resonance, in the wings of the line, or off a dust grain. 
For wing scattering the phase functi on depends on Bin as 
P(kout,eiu|wing) = |[l-(koufeiu)2] (JRvbicki fc Loebll 19991 : 



^1 



iDiikstra fc Loebl[2008. ) 

The total flux that the photon then contributes to the 
relevant pixel on the image - in this case at {xs, ys), is given^^l 

The standard factor of 47r is sometimes missing from the de- 



nominator in the literature (e.g. in lTasitsiomil2006l : lLaursen et al.l 
l200ty . because of the normalization of the phase- functions. Our 
phase functions are normalized as Ja!f2P(kiu,kout, ei^) = 47r. 



by S' = ^^^ , X P, where C — Ltot/Nj. Here, Ltot denotes 
the total Lya luminosity of the source, and A^.^ denotes the 
total number of photons used in the Monte- Carlo run (also 
see Tasjtsiomi 2006; Laursen fc Sommer-Larscn 2007 ). 



Following lRvbicki fc Loebl ( 19991 ) the linear polarization 
■p at a given location is determined by the polarized fluxes 
5*1 and Sr as 



V: 



Si — Si 
Sr + Sl' 



(A4) 



The total contribution of the photon to these polarized 
fluxes is given by 



Si= Sx [3(^)(l-cos2x) + i(l-5(M))], 
Si= Sx [5(/i)cos2x+i(l-5(M))], 



(A5) 



when a Lya photon is scattered by a hydrogen atom. 
In this expression, g(/i) = 1 for wing scattering, and 
5(m) = -.-.'"/t^ a for resonant scattering (JRvbicki fc Loebl 



ii/.i-^"^ 



1 19991 : iDiikstra fc Loebl l200d) . Furthermore, x denotes the 
angle between the photons polarization vector after scat- 
tering, denoted with Sout, and the the vector Xsca t pro- 
jected onto the x — y plane (as in iRvbicki fc Loebl Il999l : 



iDiikstra fc Loebll2008l ). In the case of wing scattering, Bout 
is obtained by finding the normalized projection of the 
old polarization vector ein onto the plane normal to kout 
iJRvbicki fc Loebl 1 19991 ). In the case of resonant scattering, 
we generate a random unit vector perpendicular to kout. 
Ea lA4l shows for example that scattering by 90° (i.e. fj, = Q) 
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Figure A4. A more detailed view of the surface brightness dis- 
tribution in the uppermost clump in the averaged image (shown 
previously in Fig |A3ll is shown. The inset shows a slice through 
this distribution along the axis x = 0. The black histogram shows 
the normalized (to the maximum) surface brightness as a func- 
tion of position y, which matches our analytic solution {red solid 
line) very well. This demonstrates that our surface brightness al- 
gorithm works well. 



results in "P = 1.0 for wing scattering, and V — 3/11 for 
resonant scattering. When a Lya photon scatters off a dust 
grain, we simply set Si — S^ — S/2. 

Figure IA3I shows the six images that we created from 
six directions (along the ±a;, y, and i-directions) , for a test 
problem in which tq = 0.1, and «c = (see above). Be- 
cause of our adopted geometry, the six images look identical 
within the noise as a result of the finite number of photons in 
the Monte-carlo run. The lower right panel shows the image 
after taking the average of all six. Figure IA4I shows a close- 
up of the upper most clump. This figure demonstrates that 
our images contain no flux where there is not supposed to 
be any. Furthermore, the surface brightness profile of indi- 
vidual clumps is in good agreement with analytic estimates: 
the inset shows a slice through the surface brightness map 
at a; = and plots the surface brightness -normalized to 
the maximum surface brightness in the clump- as a func- 
tion of y. The histogram shows the result obtained from our 
Monte-Carlo code, while the red solid line shows our analytic 
estimate, which we compute as follows: 

The total flux that we expect to detect is S{y) oc 
I-i'(y)d-sr{y,s)f{s,y). Here, f{y,s) ex (y^ + s^)"^ denotes 
the incoming flux at position {y,s), where s denotes the po- 
sition along the line of sight. This flux intersects the line of 
sight at an angle that increases with s, and the probabil- 
ity that the photon is scattered scales as r ex sjy^r?' jy- 
We can therefore write S{y) oc i /!',^^) ds (y^ -^ s^^-Va. 
At a given y, we know we will exit from the clump at 
±l(%j) = {R^-{d-yfY''^. The resulting S{y) -normalized to 
its maximum- is overplotted as the red solid line. The agree- 
ment between our analytic and Monte-carlo calculations is 
excellent which further confirms that or surface brightness 
algorithm is working well. 

We also found that the linear polarization lies in the 
range P = 23 — 27% when measured across the 4 clumps 
that are not at (a;, y) = (0, 0) in the averaged image. This 
is very close to the maximum linear polarization, Vu 



u 



analytic 

hist = MC 




Figure A5. Spectrum of photons emerging from the six outflow- 
ing, optically thick clumps. Photons were emitted in the centers of 
all clumps. The red solid line shows the analytic solution (Eq IATI I. 
while the histogram shows the solution that we obtain from our 
Monte-Carlo simulation. This plot illustrates that our code treats 
scattering in the wing accurately. 



Yi, that is expected from resonantly scattered Lya radi- 
ation. The polarization of radiation coming from the cen- 
tral clumps is consistent with in the center and rises to 
V ~ 1% on the edges, which is again consistent with analyt- 
ical expectations, which yield V — -^-^Z^ a iJDiikstra fc Loebl 
200q). For scattering in the central clumps we have sin 9 « 
tan^ « y/d, and the maximum polarization for the central 
clump is « ^ X T^ « 0.9%, which is in close agreement with 
our Monte-Carlo calculations. 



A2 Extremely Opaque Clumps with Embedded 
Lya Sources 

The previous sections showed that our code works well in 
the optically thin regime. These tests were important as 
they demonstrates the accuracy of our code with respect 
to resonant scattering, inter and intraclump propagation, 
the surface brightness and polarization algorithms. In this 
section we briefly describe one test of scattering in an ex- 
tremely opaque medium. This tests scattering in the wings 
of the line in greater detail. 

We insert photons at line center in the centers of all 
clumps. The line center optical depth through the clumps is 
enhanced to to = W^ ■ The spectrum of photo ns emerging 
from individual clumps is known analy tically (jHarringtonl 
1 19731 : lNeufeld|[l990l : iDiikstra et alll2006h and is given by 



"sph y-^ ) 



a„ToV24 ^ ^ p^jgjj 



(A6) 



This frequency x is measured in the frame of the clump. 
When a clump is outflowing, then the proper Doppler boost 
should be applied. Under the -reasonable- assumption that a 
negligible fraction of the photons scatters in a second clump, 
then the total spectrum of photons emerging from the six 
outflowing clumps 



max, res 



2A2; 



J(^) = TTT- / dx' Joph{x') 



(A7) 



X — /S.X 
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, where Aa; = «c/«th- The red solid line shows Eg IA7I for 
To = 10^, and Vc: = 100 km s^^, while the histogram shows 
the spectrum of photons that escapes from out Monte-carlo 
simulation. This plot illustrates that our code treats scatter- 
ing in the wing accurately. Further tests of our code regard- 
ing scattering i n extr emel y opaque m edia were presented in 
iDiikstra et 'sA\ l|2006h and lDiikstra~fc Loeb (2008 ). 

A3 Dust Scattering & Absorption 

Analytic expressions for the fraction of photons that escape 
from unifo rm slabs (inf i nite p lane parallel media) h ave been 
derive d by iHarringtonl l| 19731 ) and iNeufeldl (Il990l ). iNeufeldl 
l|l990 ') has provided an approximate expression for the es- 
cape fraction of Lya photons for the case in which pho- 
tons are emitted at line center in the midplane of an 'ex- 
tremely opaque' slab, where 'extremely opaque' quantita- 
tively means a^ro ^ 10^. This expressioio is 



1 F 



/esc = |cosh(3.46(a„To)'/^(l 



A)r,Y'Y 



(A8) 



where Td is the total (absorption + scattering) optical depth 
in dust from the midplane to the edge of the slab, and 
where A denotes the albedo. Figure IA6I shows that the es- 
cape fraction that we derive from our Monte-Carlo code 
agrees very well with this analytic result for zero and non- 
zero diK^_aJhedos_(as_was_aJsoJound__by_ot^ 
e.g. iLaursen fc Somm er-LarsenI 120071 : 1 Forero- Romero et al.l 
l201ll : lYaiirnret al...2012al v 



APPENDIX B: DERIVATION OF ANALYTIC 
EXPRESSIONS FOR THE SURFACE 
BRIGHTNESS & POLARIZATION PROFILES 

The surface brightness (per unit area) of scattered radiation 
at impact parameter h ± A6/2 is given by 



X /esc (a;", 6, s)sin(s, 6, x)pscat (s, 6, x)Ab. 



(Bl) 



The factor 27rfe ds Sin(s, 6, x) in the first line denotes the to- 
tal flux incident on a ring of radius h at frequency x at 
line of sight coordinate s. The probability pscat(s,fe, a:)A& 
denotes the fraction of this flux that is scattered towards 
the observer. The probability pscat(&, s,x)Ab = /c(r)A6-^ x 
(1 — exp[— rciump(2;', s, b)\), where the first term denotes the 
probability that the photon hits a clump over the path of 
length Ar = A&^, and the second term denotes the probabil- 
ity that the photon gets scattered by this particular clump. 
The incoming flux at location (6, s) and frequency x, as well 
as the expression for exp[— rciump(2:', s, 6)], is given in the 
main paper. The term f csc{x" ,b,s) denotes the fraction of 



^^ Tlie numerical factor given in INeufeldl ill 9901 ) also applies to 
our calculations, despite the different definition of tq . This is b e- 
cause the original approximate form derived by INeufeldl l l 19901 ) is 
given by/eac = l/coshCKg^' ), where Yq = [3fl4>{xs)]^/^To. Here, 
/3 = (1 — A)t^/to, and Xa = 0.525(aro)^ ". If we properly rescale 
To.neuteld ^ \/^To,us, and </'(a;)neufeld ^ </'(a;)us/\/7r, then we get 
back the original equation. 
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Figure A6. The escape fraction, /esc, of Lyo photons from a uni- 
form slab as a function of dust absorption opacity, Ta = (1 — yl)rd, 
from the midplane to the edge of the slab. Here, r^ denotes the 
total dust opacity, and A denotes the assumed albed o of the 
dust g rains. The solid line shows the analytic solution of lNeufeldl 
|l990|). The black filled circles {red filled triangles) denote the 
escape fraction obtained from our Monte-Carlo code assuming 
A = {A = 0.32). Destruction of Lyo photons is captured accu- 
rately in both cases of zero and non-zero Albedo. 




Figure Bl. This figure depicts our adopted geometry and coordi- 
nates for the analytic calculation of the surface brightness profile 
of the scattered radiation. 



the photons scattered at location (6, s) and frequency x, that 
are observed. Finally, in the main paper we expressed impact 
parameter in proper kpc, and we added a factor (kpc/asec)^ 
to convert the surface brightness proflle into the proper units 
of erg s~^ cm~^ arcsec"^. 

Equation IB4I is relevant for the total observed flux. We 
can derive the expressions for polarized Huxes Si(b) and Sr{b) 
from the scattering matrix that des cribes scattering in the 
wing of the line, which is given by (jDiikstra fc Loebll2008l . 
and references therein) 





^"2V 1 



M 



(B2) 



where the scattering matrix is defined as 
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n 

I' 



= R 



(B3) 



Here, /;' and I^ denote the components of the scattered in- 
tensity parallel and perpendicular to the plane of scattering, 
respectively. The incoming flux is directly coming directly 
from the source and is therefore likely unpolarized, in which 
case we have Ii — Ir — 4/. We can therefore write 



S,{b) 
Sr{b) 



X 



1 



4 2TvbAb 



ds 27r6 / dx X 



Xfesc{x",b,s)Sin{s,b,x)psca.t{s,b,x)AbX < 



(B4) 



APPENDIX C: DIRECTIONAL-DEPENDENT 
REDISTRIBUTION 

The frequency redistribution function, often denoted with 
-R(a;out|a;in), denotes the probability density function for the 
photon frequency a;out after scattering, given it had a fre- 
quency Sin prior to scattering. This function is an impor- 
tant quantity in the Lya radiative transfer process, and 
analytic expressions have been known for decades (see e.g. 
|Leg[l97J, and references therein). The frequency redistri- 
bution function averages over all possible scattering angels, 
/J.. However, the redistribution function varies strongly with 
yu: this is most evident when considering the case p = 1. 
Here, energy conservation implies that the photon frequency 
before and after scattering must be identical, and hence 
R{xout\xin, fi = 1) 7^ 7?(2;out|3::in). In this section we present 
a complete derivation of R{x out\x in, fJ-)- 

The photon frequencies before and after scattering are 
relateco through the angle at which the photon is scattered, 
and the total 3-D v elocity of the hydroge n atom that scatters 
the photon as (e.g. iDiikstra et al.ll2006r ) 



kin 



V ■ ko 



"th 



^th 



■+3(^-1)- 



-0{vl/c^) 



(CI) 



where g = hAuc/{2kBTc) = 2.6x 10~''(rc/10* K)'^^^ is the 
fractio nal amount of energy that is transferred per scattering 
event l|Fieldl 1 19591). Throu ghout this calculation we safely 
ignore recoil iJAdamdl 19711 ). 

For simplicity, but without loss of generality, we can de- 
fine a coordinate system such that kin = (1, 0, 0), and kout = 
(/i, ^1 — fi^, 0), i.e. the photon wav e vector s lie entirely in 
the x-y plane. Following lAhn et al.l (|2000|) we decompose 
the atom's velocity into components parallel (nn) and or- 
thogonal {vy and Vz) to kin, and we have v — {v^^,Vy,Vz). 



Eq ICll can then be recast as 






(C2) 



Xin — U + Ufl + W V 1 — H'^, 

where we have introduced the dimensionless velocity 

^* This assumes coherence in the frame of the atom, which is 
relevant at the densities and temperatures of interest (see e.g. 
Dijkstra et al. 2006 for a more detailed discussion). 



parameters u = D||/vth and w = Vy/vth- Note that the value 
of Vz is irrelevant in this equation. 

The velocities u and w are unrelated, and the most gen- 
eral way of writing the directional redistribution function is 

du / dw (C3) 

-OO J — CXD 

R{Xout\lJ., a;in, U, w)P{u\^i, Xin)P{w\fl, a;in), 

where A/" denotes the normalization constant. The integral 
over w can be eliminated by utilising Eq IC2I That is, we 
replace R{xout\fi,Xin,u,w) = 5_d(/[wu]) = Soiw - Wu)/-^, 
in which f[wu] = a;out — Sin + u — ujj. — Wu\/^ — A*^. We find 



R{Xout\^J.,Xin) =Af duP{u\^,Xin)P{Wu\fJ.,Xin), (C4) 



where the factor 1/^ = l/\/l — fi^ is absorbed by the 
normalization constant J\f. 

The conditional absorption probabilities for both w 
and u cannot depend on the subsequent emission direction, 
and therefore P{u\^, Xin) — P{u\xin) and P{wu\^J.,Xin) — 
P{wu\xin). Furthermore, the conditional PDF for w does 
not depend on lin either. This is because w denotes the 
normalized velocity of the scattering atom in a direction 
perpendicular to kin, and the frequency that the atoms 
'sees' does not depend on w. The absorption probabil- 
ity can therefore not depend on w, and P{'Wu\xin) ~ 

P{wu) = \/ 2TTk'' T 6xp(— w^), where we assumed a Maxwell- 
Boltznrann distribution for the atoms' velocities. 

The expression for P(u\x i n) ca n be obtained from 
Bayes theorem (see e.g. iLed |l97J), which states that 

P{u\xin) = P{u,Xin)/P{xin) = P{Xin\u)P{u) / P(xin) , m 

which P{xin\u) denotes the absorption probability for a sin- 
gle atom that has a speed u, and P(a;in) oc a{xin). There- 
fore, P{u\xin) OC P(a;in|u)P(u)/a(3;in), where P{xin\u) — 

^ [^„fe„-^Kh/c]^+Ai^/4 (e.g. i Rybicki fc Lightmanl Eili). 
If we substitute this into Eo IC4I and absorb all factors that 
can be pulled out of the integral into the normalization con- 
stant A/", then we get 



R{Xout\fJ., Xin) = M" du 



exp(—u 



xexp 



(Sin — u)2 -f ag 
fAx + u{^- 1U2 



(C5) 



r /Ax+u 



, where we introduced Aa; = Xi^ — x^y^x.- The normalization 
constant can be computed analytically, and we have 



- T^-^..., V - r- 



du, ""P^T!'^ . (C6) 



(a;in — u)2 + ( 



\- 



fAx + u{fl- 1)\2 



xexp I - (^ /)—^ ' ' I ^" '^' "^ ^• 



For /i = 1 we have a;out = a;in fEg IC2[) . For fi — ~1 we 
have a^out ~ xin — 2u, and we have _R(a;out|/i = — l,a;in) = 
|P(uc|a;in), in which Uc ~ (sin — a:ont)/2. We have verified 
that these analytic expressions are in excellent agreement 
with results obtained from Monte-carlo calculations. 
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